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SUMMARY 

A  detailed  description  is  given  of  a  technique  for  converting  a  contour  map  into 
an  equispaced  grid  of  points.  A  full  description  is  given  of  program  MAPGRID  which 
converts  digitized  contour  data  into  such  a  grid  of  points. 
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I  INTRODUCTION 

I  .  I  Ob j  ect ive 

This  Report  describes  a  technique  for  converting  a  contour  map  into  an  equispaced 
grid  of  points.  The  purpose  of  the  grid  representation  is  to  enable  the  map  data  to  be 
handled  more  easily  in  a  computer  and,  as  such,  it  is  to  form  the  base  on  which  all 
operations  on  the  map  data  are  to  be  performed.  It  is  assumed  that  the  contours  repre¬ 
sent  some  physical  quantity  and  thus  that  the  variable  represented  by  the  contour  height 
is  a  single-valued  function  of  the  map  co-ordinates.  It  is  also  assumed  that  the 
variable  is  smooth  and  continuous. 

The  technique  described  here  has  applications  in  image  processing  where,  for 
instance,  it  may  be  necessary  to  compare  (by  ratioing)  an  image  and  a  contour  map.  A 
grid  representation  is  also  an  essential  prerequisite  for  degrading  the  spatial  resolu¬ 
tion  of  a  contour  map  by  convolution,  thus  enabling  contour  maps  having  different 
resolutions  to  be  compared  at  a  common  resolution.  These  and  other  applications  of  this 
technique  are  to  be  described  in  a  subsequent  report. 

Included  in  this  Report  are  discussions  of  the  problems  associated  with  the  digiti 
nation  of  a  contour  map  and  those  of  interpolation  in  general.  The  method  used,  the 
interpolation  algorithms  and  the  reasons  for  their  use  are  described  in  detail.  Also 
described  are  two  operations  which  may  be  required  to  be  performed  on  a  grid;  those  of 
expanding  (or  contracting)  the  number  of  data  points  in  the  grid,  and  the  construction 
of  cross-sections.  Listings  of  all  programs,  which  are  written  in  computer  independent 
ANSI  FORTRAN  V,  are  given. 

I . 2  Outline  of  main  program 

Program  MAPGRID  transforms  digitized  contour  data  into  a  rectangular  grid  of 
points.  The  program  takes  two  orthogonal  sets  of  cuts  across  the  contour  map  (see  Fig  1 
the  points  of  intersection  between  these  grid  lines  define  the  points  at  which  interpo¬ 
lated  data  is  required.  For  each  grid  line  the  points  of  intersection  of  the  line  with 
contours  are  calculated  and  then  interpolation  is  carried  out  along  the  grid  line  to 
obtain  values  at  the  grid  points  lying  on  the  line.  Also  produced  are  weighting  values 
reflecting  the  confidence  in  the  interpolated  values.  This  procedure  results  in  two 
grids  of  points  being  ohtained  corresponding  to  the  two  sets  of  cuts,  together  with  two 
grids  containing  the  associated  weighting  values.  These  weighting  values  are  used  to 
combine  the  grids  to  produce  the  final  grid. 

It  is  necessary  to  take  two  sets  of  cuts  since  it  may  he  that  some  contours  run 
parallel  to  one  of  the  sets  of  cuts  and  thus  would  be  ignored  by  that  set  unless  they 
fell  exactly  on  grid  lines.  Contours  missed  by  one  set  of  cuts  are  therefore  incor¬ 
porated  in  the  orthogonal  set.  In  this  way  the  orientation  of  the  cuts  with  respect  to 
the  contour  map  is  not  important,  whilst  it  may  be  crucial  if  onlv  one  set  of  cuts  were 
to  be  taken. 
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2  DIGITIZATION 

2. I  Digitization  of  contours 

In  general,  the  investigator  will  not  have  the  contour  information  already  digi¬ 
tized,  but  rather  will  only  have  a  paper  copy  of  the  imd,  thus  necessitating  the  digi¬ 
tization  of  the  contours.  There  are  several  reasons  why  the  contours  should  be  digitized 
rather  than  digitizing  along  the  grid  lines  that  the  computer  program  is  going  to  use: 

(i)  It  is  very  difficult  to  determine  accurately  the  point  of  intersection  of  a 
grid  line  with  a  contour  if  the  two  cross  at  a  very  oblique  angle. 

(ii)  Digitization  around  a  contour  is  often  much  easier  if  a  manual  digitizing 
table  is  to  be  used,  since  the  exact  position  of  digitized  points  along  a  contour  is  not 
important,  whereas  the  exact  point  of  intersection  between  a  grid  line  and  a  contour  is 
much  more  important. 

(iii)  If  in  the  future  it  is  required  to  generate  grids  at  other  larger  spacings, 
then  it  is  not  necessary  to  redigitize  the  map. 

(iv)  It  is  easier  to  check  for  any  errors  in  the  digitized  data  if  the  digitiza¬ 
tion  is  performed  around  contours,  because  the  contours  can  easily  be  plotted  and  com¬ 
pared  with  the  original.  With  the  alternative  method  a  complex  contour  threading 
operation  would  have  to  be  performed  before  the  contours  could  be  plotted,  which  would 
in  itself  give  rise  to  some  differences  with  the  original  contours. 

It  is  necessary  for  the  investigator  to  decide  the  interval  at  which  points  on  a 
contour  should  be  digitized.  The  evaluation  of  this  interval  involves  an  estimate  being 
made  of  the  interval  along  a  contour  that  can  be  represented  by  a  straight  line.  This 
interval  is  dependent  upon  the  grid  spacing  and  the  curvature  of  the  contour  in  question; 
the  smaller  the  grid  spacing  or  the  greater  the  curvature  of  the  contour,  the  smaller 
this  interval  must  be  in  order  to  represent  the  contour  to  the  required  accuracy.  More 
precisely,  if  d  is  the  digitization  interval,  r  is  the  radius  of  curvature  of  the 
contour  (which  generally  is  a  function  of  position  along  the  contour)  and  s  is  the 
maximum  permissible  deviation  of  a  straight  line,  joining  two  points  on  the  contour, 
from  the  contour,  then  d  %  /8r s  ,  where  d  is  such  that  the  angle  the  contour  turns 
through  is  small.  The  author  has  found  that  for  his  work  it  is  adequate  that  s  be 
half  the  grid  spacing,  but  this  may  not  be  acceptable  for  all  applications. 

If  the  digitization  interval  is  greater  than  d  then  the  contour  is  under¬ 
digitized,  and  if  it  is  less  it  is  over-digitized.  The  only  consequence  of  the  latter 
is  that  the  number  of  data  points  describing  the  contour  is  greater  than  necessary.  The 
contour  fitting  algorithm  (described  in  section  3.3.3)  is  considerably  better  than 
linear  interpolation  and  is  adopted  to  make  a  'best'  estimate,  which  is  necessary  with 
under-digitized  data. 

The  digitization  of  contours  can  take  either  of  the  following  forms  (amongst 
others) : 

(i)  Constant  separation  between  data  points. 
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(ii)  Variable  density  system  with  the  interval  depending  upon  the  curvature  of 
the  contour.  This  reduces  the  number  of  points  required  to  describe  a  contour. 

If  it  is  decided  to  opt  for  a  constant  separation  between  data  points  then  the  map 
will  be  correctly  digitized  if  the  evaluation  of  the  digitization  interval  is  based  upon 
the  smallest  radius  of  curvature  of  any  part  of  any  contour.  The  variable  density 
method  requires  a  continual  evaluation  of  the  curvature  of  the  contour  and  is  best 
suited  to  automated  digitization.  If  the  digitization  is  being  performed  manually  using 
this  method  then  the  investigator  must  make  a  subjective  estimate  of  the  required 
digitization  interval. 

2. 2  Grid  spacing 

The  main  information  content  of  a  contour  map  is  contained  in  the  relative  posi¬ 
tions  of  the  different  contours  and  the  large  scale  features  shown  by  them,  rather  than 
in  any  small  scale  detail  that  may  be  shown  by  individual  contours.  Thus  to  retain  the 
main  information  content  when  transforming  a  contour  map  into  a  grid,  the  grid  spacing 
must  be  half  the  value  of  the  minimum  separation  between  any  two  contours  or  of  large 
scale  loops  of  the  same  contour.  Any  detail  on  a  scale  smaller  than  this  will  not  be 
contained  in  the  grid.  More  generally,  the  grid  spacing  should  be  at  least  half  the  size 
of  the  smallest  feature  of  interest. 

2. 3  Author's  note 

The  above  discussion  is  intended  to  illustrate  some  of  the  basic  problems  involved 
in  digitizing  contours.  It  is  not  intended  to  be  an  exhaustive  or  highly  rigorous 
treatment  of  the  subject. 

3  INTERPOLATION  TECHNIQUES 

3. 1  The  problem  of  interpolation 

It  is  often  desired  to  estimate  the  value  of  a  function,  which  is  sampled  at  cer¬ 
tain  discrete  points,  at  some  intermediate  points.  This  process  is  known  as  interpola¬ 
tion.  Although  any  set  of  data  points  can  be  interpolated  by  an  infinite  number  of 
different  functions,  one  normally  requires  that  interpolated  values  are  reasonable.  What 
is  believed  to  be  reasonable  depends  upon  the  laws  and  processes  underlying  the  data 
values.  It  must  be  realised  that  unless  a  precise  mathematical  function  expressing  those 
laws  can  be  formulated,  then  no  method  of  interpolation  can  accurately  reproduce  the 
missing  data  points.  In  many  instances  it  may  be  that  whilst  the  nature  of  the  processes 
are  understood  they  are  not  expressed  in  mathematical  form,  eg  surface  topography  con¬ 
sidered  as  an  outcome  of  geological,  climatic  and  cultural  factors.  In  such  instances  it 
is  not  possible  to  distinguish  analytically  and  quantitatively  between  the  desirability 
of  different  algorithms.  Instead,  the  choice  has  to  rely  on  the  judgement  of  the  inves¬ 
tigator  himself,  ?'<;  some  rationalisation  of  a  comparison  between  an  achieved  interpola¬ 
tion  value  and  a  subjective  impression  of  what  the  value  would  have  been  at  that  point. 

In  cases  where  the  nature  of  the  underlying  processes  and  laws  are  not  known  the  normal 
criteria  of  smoothness  and  simplicity  have  to  be  adopted.  These  are  conditions  which  are 
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iituitivei;  desirable  in  the  absence  of  detailed  information  as  to  the  nature  of  the  fit 
required. 

Interpolation  functions  can  be  of  two  types:  local  or  global.  A  local  function 
uses  only  data  points  near  to  the  point  at  which  an  interpolated  value  is  required, 
eg  linear  interpolation.  A  global  function,  however,  utilizes  all  of  the  data  points  in 
the  data  set,  eg  cubic  spline  interpolation.  With  this  function  distant  data  points  make 
a  contribution  to  the  interpolated  value.  Also,  if  the  data  set  contains  a  large  number 
of  points  then  this  type  of  function  generally  involves  a  large  amount  of  computation. 

In  cases  where  there  is  a  mathematical  function  expressing  the  variable  it  may  be  that  a 
global  function  is  the  most  appropriate,  but,  in  instances  where  the  only  requirements 
are  for  smoothness  and  simplicity,  then  a  local  function  can  be  adopted.  Indeed,  the 
adoption  of  such  a  requirement  is  desirable  if  there  is  no  knowledge  of  the  underlying 
processes;  the  adoption  of  a  global  function,  whereby  distant  data  points  can  affect  the 
interpolated  values,  unnecessarily  implies  more  correlation  between  data  points  than  is 
warranted.  Obviously  if  it  is  known  that  there  is  no  correlation  between  near  and  dis¬ 
tant  data  points  then  a  local  function  must  be  adopted. 

The  suitability  of  a  particular  function  is  also  dependent  upon  the  density  of  the 
data  points.  For  very  closely  spaced  data  the  nature  of  the  interpolation  function  is 
largely  unimportant  and  it  may  be  that  linear  interpolation  is  entirely  appropriate.  In 
general,  the  more  widely  spaced  the  data  the  fewer  functions  will  prove  acceptable.  A 
further  (practical)  consideration  is  the  amount  of  computation  involved  with  different 
functions  if  large  numbers  of  data  points  and  interpolated  values  are  required.  Thus  a 
compromise  may  have  to  be  made  between  the  suitability  of  a  particular  algorithm  and 
maintaining  the  amount  of  computation  at  a  manageable  level. 

In  choosing  an  interpolation  function  it  should  be  realised  that  in  general  no 
method  of  interpolation  is  precisely  accurate  except  at  the  points  through  which  the 
function  has  been  fitted.  If  little  is  known  about  the  exact  properties  of  the  variable 
involved,  then  the  assessment  of  the  suitability  of  any  interpolation  function  is  largely 
left  to  the  investigator's  intuition.  Thus  no  general  rule  as  to  the  applicability  of  a 
particular  function  can  be  given,  so  each  problem  must  be  considered  individually. 

3.2  Interpolation  requirements 

There  are  two  different  interpolation  procedures  that  have  to  be  carried  out: 

(i)  Interpolation  along  a  contour. 

(ii)  Interpolation  along  a  grid  line  using  contour  cut  data. 

As  mentioned  earlier,  the  choice  of  an  interpolation  function  is  dependent  upon  the 
density  of  the  data.  For  interpolation  around  a  contour,  linear  interpolation  will  be 
appropriate  if  the  spacing  between  the  data  points  is  sufficiently  small;  however,  if 
this  is  not  the  case,  a  more  appropriate  function  is  required.  This  function  should 
possess  the  following  properties: 
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(i)  Continuous. 

(ii)  Continuous  first  derivative. 

(iii)  Local. 

(iv)  Independent  of  the  co-ordinate  system. 

The  first  two  requirements  are  those  normally  adopted  where  no  knowledge  of  the 
underlying  laws  and  processes  is  assumed  and  they  are  appropriate  here,  because  contours 
are  expected  to  be  smooth  and  continuous.  Since  it  is  to  be  taken  there  is  no  a  priori 
knowledge  of  the  properties  of  the  variable  involved,  it  is  required  that  the  interpola¬ 
tion  algorithm  be  local  rather  than  global,  ie  dependent  only  upon  the  data  points  near 
to  the  points  at  which  interpolated  values  are  required.  This  also  means  that  the 
interpolation  function  is  going  to  be  reasonably  efficient  for  implementation  on  a  com¬ 
puter.  The  fourth  requirement  is  introduced  because,  for  contours,  the  co-ordinate 
system  (generally)  does  not  have  any  significance  for  the  variable  involved,  thus  the 
shape  of  an  interpolated  curve  should  be  independent  of  the  orientation  of  the  axes  of 
the  map  co-ordinate  system.  The  choice  of  a  suitable  function  for  contour  data  is 
further  complicated  in  that,  in  cartesian  co-ordinates,  infinite  gradients  (dy/dx  =  «°) 
can  be  encountered,  and  the  function  may  have  to  be  multi-valued  under  certain 
circumstances. 


When  interpolating  intermediate  values  along  grid  lines,  infinite  gradients  and 
multi-valuedness  will  not  be  encountered  and  there  is  no  need  for  the  function  to  be 
independent  of  the  orientation  of  the  co-ordinate  system  because  the  choice  of  axis  for 
the  contour  'height'  is  not  arbitrary;  no  other  choice  of  axis  could  reasonably  be 
adopted  for  physically  realisable  data.  Thus  for  interpolation  along  a  grid  line  a 
function  is  required  that  is  continuous,  with  a  continuous  first  derivative  and  is 
local.  In  addition,  the  function  must  be  constrained  at  turning  points  (ie  maxima  and 
minima),  so  that  interpolated  values  do  not  reach  a  neighbouring  contour  level. 

With  these  considerations  in  mind,  it  was  decided  to  adopt  an  average  tangents 
form  of  interpolation. 

3.3  Average  tangents  interpolation 


3.3.1  The  basic  algorithm 

The  method,  which  is  applicable  to  any  data  set  for  which  y  is  a  single-valued 
function  of  x  ,  is  illustrated  in  Fig  2.  T^  is  the  connection  vector  between  adjacent 
data  points  and  the  tangent  of  the  angle  between  the  positive  x  axis  and  T^  is  given 
by: 


tan(e  . ) 


i  =  I,  2, 


I 


where  n  is  the  number  of  data  points  in  the  data  set. 
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An  average  is  now  taken  of  pairs  of  consecutive  tangents  as  follows: 

tan (6 .  , )  +  tan (0 . ) 

A.  =  - ~~~~2 - ~  i  =  2,  3,  . .  n  -  1 

For  the  first  and  last  intervals: 


A^  =  tan  (m ^ ) 


A  =  tan(y  ) 
n  n- 1 


The  average  tangent  A^  represents  the  mean  gradient  at  the  ith  data  point.  Since 
an  interpolation  function  is  required  that  is  continuous  with  a  continuous  first  deriva¬ 
tive  the  function  can  he  represented  by  a  third  order  polynomial  between  points  i  and 
i  ♦  1  »  ic 

3  2 

y  =  a(x  -  x.)  +  b(x  -  x.)  +  c(x  -  x.)  +  d  . 

The  coefficients  of  the  polynomial  can  be  found  by  choosing  the  first  derivative 
of  the  polynomial  at  the  points  i  and  i  +  1  as  being  the  values  of  the  average  tan¬ 
gent  at  these  points,  whence: 

Ai(Vi  ■  xi)  +  Ai+I(xi+1 "  V  '  2(yi+i  ■  yi> 
a  _  - 

(xi+l  -  X.) 


b  = 


3(yi+l  -  yj)  -  2Ai(xi*l  -  xi>  -  Aj+l(Vl  -  xj> 

(xi+i  -  x.)2 


c  =  A. 

l 

d  =  yi  • 


With  this  method  a  different  third  order  polynomial  is  fitted  to  each  interval. 


3.3.2  Constraining  the  interpolation  function  near  turning  points 

The  interpolation  function  for  contour  cut  data  has  to  be  modified  so  that  inter¬ 
polated  values  remain  within  the  interval  defined  by  the  two  points  between  which  the 
function  is  valid,  or,  if  this  interval  is  zero,  do  not  reach  a  particular  value. 

Fig  3aSb  (bold  lines)  illustrates  cases  where  such  modification  is  required.  The  x  and 
y  values  at  the  turning  points  (denoted  by  subscript  t  )  are  obtained  from  the 
equations : 


and 


2 

0  =  3ax  +  2bx  +  c 

t  t 

3  2 

y  =  ax  +  bx  +  cx  +  d 
7t  t  t  t 


If  both  values  of  xf  lie  between  x.  and  x.+]  then  the  situation  is  as 
illustrated  in  Fig  3a,  and  if  only  one  value  of  xf  satisfies  this  criterion  then  the 
situation  is  as  depicted  in  Fig  3b.  The  former  case  will  be  discussed  first. 
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If  both  values  of  y  lie  between  y. 


and 


y^+^  the  interpolation  function  is 


t  '  1 

not  modified;  otherwise  the  range  is  divided  into  three,  defined  by  the  points  given 
below,  and  new  interpolation  functions  fitted. 


X  =  X. 

1 

A 

=  A. 

l 

y 

=  yi 

X  =  xtl 

A 

=  0 

y 

=  ytl 

if 

ytl 

< 

yi+l 

y 

=  y.  +  a(y.+l  -  y.) 

if 

yt| 

> 

yi+l 

X  =  Xt2 

A 

=  0 

y 

=  yt2 

if 

yt2 

> 

yi 

y 

=  yi+l  "  “(yi+l  '  yi} 

if 

yt2 

* 

yi 

X  =  X  .  , 

1  +  1 

A 

=  A.  , 
i+l 

y 

=  yi+l 

(a  is  a  positive  constant  less  than  1  such  that 


y  <  for  the  precision  with 


which  interpolated  values  are  required.  Generally  a  =  0.99  is  a  suitab’e  choice.) 


If  the  new  interpolation  function  does  not  remain  between  y.  and  y.  , 

'l  l+l 


(as 

values  near  the 


illustrated  in  Fig  3a  by  the  dashed  line)  then  the  interpolated  y 

turning  point  that  lie  beyond  either  y.  +  a(y.  ,  -  y.)  or  y.  ,  -  a(y.  ,  -  y.)  , 

’  l+l  l  'l+l  l+l  l 

whichever  is  appropriate,  are  set  to  that  value.  This  situation  will  arise  if 


■■  ■feR)  ■ 


a  condition  that  for  most  types  of  data  should  rarely  occur.  The 


procedure  just  described  gives  rise  to  discontinuities  in  the  gradient,  however  in 
practice  since  interpolated  values  are  required  at  discrete  points  this  somewhat 
alleviates  the  problem. 


If  there  is  only  one  value  of  x^  in  the  range  x^ 


to 


x.+^  it  may  be  required 


y  (ie  the 
o 


that  interpolated  values  are  not  to  exceed  a  particular  value  of  y  ,  say 
next  contour  level).  If  the  interpolation  function  does  not  reach  y  then  the  func¬ 
tion  remains  unchanged,  otherwise  two  new  functions  are  fitted  using  the  following 
points: 


x . 
l 


x . 
l  +  l 


A  = 


A. 

l 


A  =  A. 


i  +  1 


v .  + 

'  l 


yi  +  l 


*<y. 


y . ) 

i 


If  the  new  function  does  not  remain  between  y.  and 

l 


v . 

■  l 


i  (v 


v  .  )  (as  i I lus- 

i 


trated  in  Fig  3b  by  the  dashed  line),  then  interpolated  values  lying  beyond 

y.  +  n(y  -  y.)  are  set  to  that  value.  It  can  be  shown  that  for  contour  cut  data, 

1  O  1 


where 


y^+j  ,  the  interpolation  function  must  take  the  form  shown  in  Fig  3b 


(ie  remain  above  or  below  y.  between  x, 

l  l 


A .  ,  have 
l+l 


and  x.  ,  )  since  A.  and 
l  +  l  i 

opposite  senses  (if  one  is  positive  the  other  is  negative;  zero  can  be  regarded  as  posi¬ 
tive  or  negative  depending  upon  what  is  necessary  for  A.  and  A.  ,  to  be  of  opposite 

i  i  +  l  rl 

sense) . 


The  average  tangents  algorithm  is  not  independent  of  the  orientation  of  the 
co-ordinate  system  and  in  some  cases  the  algorithm  fails  completely,  j  when  infinite 
gradients  are  involved,  when  y  is  a  multi-valued  function  of  x  .  This  problem  can 
however  be  overcome.  Rather  than  considering  y  changing  with  respect  to  x  , 
x  and  v  are  considered  independently  as  changing  with  respect  to  the  separation  (s) 
between  adjacent  data  points. 

Def ine 


The  average  tangents  algorithm  can  now  be  applied  to  x  and  y  separately  as  a 
function  of  with  parameters  m.  being  substituted  for  tan(0^)  .  Two  equations 

are  thus  obtained: 


and 


x 


3  2 

as  +  bs  +  cs+d 
X  X  X  X 


3  2 

y  =  as  +  bs  +  cs  +  d 

y  y  y  y 

where  s  =  ^(x  - 

Thus  to  derive  the  y  value  corresponding  to  a  particular  x  value,  the  value  of 
s  corresponding  to  the  x  value  is  found  by  solving  the  first  equation  (with  the  condi¬ 
tion  that  0  «  s  «  s.  ),  and  then  substituting  the  value  for  s  in  the  second  equation 
to  obtain  y  , 

With  this  procedure  the  algorithm  never  fails  because  'infinite'  gradients  are 
never  encountered  since  s.  >  0  (s.  =  0  corresponds  to  two  adjacent  data  points  being 

coincident).  The  algorithm  is  independent  of  the  orientation  of  the  co-ordinate  system 
since  x  and  y  are  treated  in  exactly  the  same  manner  as  a  function  of  the 
co-ordinate  free  quantity  s  . 

4  PROGRAM  DESCRIPTION 


-  * 1 


Details  of  the  main  program  such  as  dimensions  of  arrays,  initialisation  of  para¬ 
meters,  and  program  units  called  are  given  in  Appendix  A.  Detailed  specifications  of 
all  subroutines  are  given  in  Appendix  B,  and  a  complete  program  listing  is  given  in 
Appendix  C.  In  the  program  listing  the  parameters  initialised  at  the  start  of  the  main 
program  are  those  for  the  example  discussed  in  this  Report. 


1 1 

4.1  Program  ptiens 

Limits  to  the  allowed  values  of  interpolated  map  data  have  to  be  set.  If  no  limits 
are  to  be  placed  on  interpolated  values  then  the  parameters  setting  the  limits  should  be 
set  so  that  they  are  well  away  from  any  of  the  contour  values  and  anticipated  interpo¬ 
lated  values.  The  parameter  FRACTN  sets  how  closely  interpolated  values  can  approach  the 
next  contour  level  at  turning  points.  It  is  recommended  that  this  be  set  at  0.99,  as  dis¬ 
cussed  in  section  3.3.2,  so  that  interpolated  values  just  fall  short  of  the  next  level. 

Maps  may  be  classified  into  two  categories  for  the  problem  being  considered  here, 
closed  or  open.  A  closed  map  is  defined  to  be  one  where  all  contours  start  and  finish 
within  the  region  being  considered  (as  in  Fig  1)  and  an  open  map  is  one  where  only  some 
contours  start  and  finish  within  the  region  (as  in  the  area  enclosed  by  the  dashed  rec¬ 
tangle  in  Fig  1).  The  operation  of  the  program  is  slightly  different  for  the  two  types 
of  map,  this  being  related  to  combining  the  values  obtained  from  the  two  sets  of  grid 
lines  and  with  the  problem  of  interpolation  outside  a  dataset.  These  are  discussed  in 
more  detail  in  sections  4.4  and  4.6  respectively.  The  two  modes  are  selected  as  follows: 

(i)  Closed  maps.  IROUTE  =1,  ZE  =  value  to  which  all  points  beyond  the  outer¬ 

most  contour  are  to  be  set. 

(ii)  Open  maps.  IROUTE  +  I.  ZE  =  value  to  which  all  points  through  which  pass 
x  and  y  direction  grid  lines  that  do  not  intersect  any  contours. 

4 . 2  Contour  data  formatting 

The  program  requires  that  the  contour  data  be  in  a  specific  form.  Each  contour 
should  have  a  header  which  should  then  be  followed  by  the  x  and  y  co-ordinates  of 
each  point,  there  being  one  point  per  line.  The  header  contains  three  parameters:  the 
contour  value,  the  number  of  points  in  the  contour  and  the  type  of  contour.  The  first 

two  par;imeters  are  self-explanatory  but  the  third  requires  further  comment. 

Contours  can  be  of  two  types,  either  closed  such  as  a  circle  or  open  such  as  a 
line.  It  is  necessary  to  know  which  of  these  types  a  contour  is  because  of  the  problems 
associated  with  interpolation  in  the  end  intervals  of  a  dataset.  The  average  tangents 
algorithm  (section  3.3.1)  sets  the  average  tangent  at  the  end  data  points  to  be  the 
gradient  between  the  two  points  forming  the  end  interval,  but  for  the  case  of  a  closed 
contour  a  better  estimate  is  possible  because  the  data  points  beyond  the  end  interval  of 
the  dataset  are  known.  Thus  for  closed  contours  (for  which  it  is  assumed  that  the  first 
and  last  points  have  the  same  co-ordinates)  the  first  two  points  are  also  added  to  the 
end  of  the  dataset  and  then  interpolation  is  only  considered  between  the  second  and 
second  to  last  points  on  the  contour.  This  procedure  cannot  of  course  be  carried  out 
with  open  contours,  yet  it  would  be  convenient  if  both  types  could  be  treated  in  the  same 
manner.  Thus  for  open  contours  an  extra  data  point  is  added  at  each  end  using  linear 
interpolation  so  that  the  average  tangents  at  what  are  now  the  second  and  second  to  last 
points  remain  the  same.  Therefore,  after  a  contour  has  been  read  in,  it  is  reformatted 
in  one  of  the  above  ways  by  subroutine  IIJENT. 
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4 . 3  Program  operation 

The  program  sLarts  by  initialising  all  the  constants  and  calculated  the  x  and  y 
values  of  the  grid  lines.  The  values  of  the  points  in  the  grid  produced  by  taking  cuts 
in  the  y  direction  are  then  calculated. 

Contours  are  read  in  and  processed  one  at  a  time.  The  contour  data  is  first 
checked  by  subroutine  CHECK  to  ensure  that  no  two  adjacent  points  have  the  same 
co-ordinates,  which  would  cause  the  contour  interpolation  algorithm  to  fail.  The 
contour  data  is  then  reformatted,  as  described  above,  by  subroutine  IDENT.  Following 
this,  subroutine  COORDS  is  called  to  calculate  the  points  at  which  the  y  direction 
grid  lines  intersect  the  contour;  the  algorithm  described  in  section  3.3.3  is  used  here. 
This  procedure  is  repeated  for  each  contour,  after  which  subroutine  HEIGHT  is  called  to 
sort  the  different  contour  values,  together  with  the  upper  and  lower  limits  to  the  map 
data,  into  ascending  order.  Each  of  the  y  direction  grid  lines  is  now  taken  in  turn. 
Subroutine  SELECT  is  called  to  obtain  all  the  points  of  intersection  of  contours  with  a 
chosen  grid  line  and  these  points  are  then  arranged  into  ascending  order  in  the  y 
direction  by  subroutine  ORDER.  Subroutine  CHECK2  is  then  called  to  ensure  that  no  two 
adjacent  data  points  have  the  same  y  value,  which  would  cause  the  subsequent  inter¬ 
polation  algorithm  to  fail.  Subroutine  AVTAN  calculates  the  values  on  the  grid  along 
the  grid  line  using  the  average  tangents  algorithm  described  in  sections  3.3.1  and 
3.3.2.  The  associated  weighting  factors  for  each  of  the  grid  points  are  obtained  by  a 
call  to  subroutine  WEIGHT. 

The  above  procedure  is  repeated  for  the  grid  lines  in  the  x  direction,  the  only 
differences  being  that  subroutine  HEIGHT  is  not  called,  and  as  each  new  line  is  gener¬ 
ated  it  is  combined  with  the  first  grid  to  produce  the  final  grid.  On  completion  of 
this  the  final  grid  is  then  written  to  a  file  with  the  data  being  written  from  the  top 
of  the  grid  downwards,  ie  as  a  series  of  x  direction  grid  lines  in  descending  values 
of  y  . 

4 . 4  Producing  the  final  grid 

Associated  with  each  interpolated  value  is  a  weighting  factor.  The  weighting 
factor  reflects  the  fact  that  an  interpolated  value  is  most  likely  to  be  reliable  where 
the  original  data  is  most  closely  spaced.  The  weighting  factor  is  given  by: 


where  the  point  at  which  an  interpolated  value  is  required  lies  between  x^  and  x-+j  • 
If  the  point  lies  outside  the  range  of  the  dataset  then  W  =  0  . 

Three  methods  of  combining  the  two  grids  were  investigated: 


(G 


G  ) 


(i) 
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(ii)  fv 


(GW  +  G  W  ) 
xx  y  y 

(W  +  W  ) 
x  y 

(G  +  G  ) 

x  y 


(W  +  W  >0) 
x  y 


(W  +  W  =0) 
x  y 


(iii)  G , 


=  G 


(G  +  G  ) 
x  y 


(W  >  W  ) 
x  y 

(W  <  W  ) 
x  y 


(W  =  W  ) 
x  y 


where  G^  and  G^  are  the  values  at  a  particular  grid  location  produced  by  taking  cuts 

parallel  to  the  x  and  y  axes  respectively,  W  and  W  are  the  corresponding 

x  y 

weighting  factors,  and  G^.  is  the  final  value. 

The  desirability  of  each  of  these  methods  was  judged  on  the  difference  between  the 
grid  obtained  from  a  series  of  contours  of  a  two  dimensional  Gaussian  and  a  grid  derived 
analytically.  The  contours  were  produced  at  intervals  of  25  between  25  and  225  and  both 
grids  were  set  to  0  beyond  the  outermost  contour.  Figs  4  and  5  show  the  two  grids  pro¬ 
duced  by  the  two  sets  of  cuts  and  Figs  6,  7  and  8  show  the  differences  between  the  pro¬ 
gram  grid  and  the  analytic  grid  for  the  three  methods;  differences  greater  than  | 9 | 
are  set  to  9. 

Fig  6  lends  weight  to  the  reasoning  for  the  weighting  factors  used  in  methods  2  and 
3,  since  the  differences  are  a  minimum  45°  away  from  the  X  and  Y  axes,  where  the 
spacing  of  the  original  data  is  the  same  for  each  set  of  grid  lines.  Fig  7  shows  a  con¬ 
siderable  improvement  in  the  agreement  but  suggests  that  considerably  greater  weighting 
should  be  given  to  the  larger  of  the  two  weighting  values.  This  is  taken  to  the  extreme 
in  the  third  method,  producing  a  further  improvement. 

It  is  the  third  method,  with  some  modification,  that  has  been  adopted  to  combine 
the  values  produced  by  the  two  sets  of  grid  lines.  The  weighting  values  are  calculated 

as  before,  except  that  if  A.  =  A.  ,  =  0  and  z.  =  z.  ,  then  W  =  -1  ,  and  if  the  grid 

l  i+l  i  l+l 

line  under  consideration  does  not  intersect  any  contour  then  W  =  -2  .  The  method  for 
combining  the  values  is  dependent  upon  the  type  of  map,  and  is  as  carried  out  according 
to  the  following  sets  of  rules: 
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The  extra  conditions  arise  from  the  following  considerations.  In  Fig  9a  the  grid 
line  yy  intersects  the  contour  with  height  2  at  points  a  through  e  .  Interpolated 
values  along  yy  between  b  and  d  are  given  as  2  but,  by  inspection,  the  variable 
should  be  less  than  2  between  b  and  c  and  greater  than  2  between  c  and  d  .  The 
correct  sense  for  the  interpolated  data  is  obtained  by  using  interpolated  data  produced 
by  the  orthogonal  set  of  grid  lines.  If  the  weighting  of  a  point  is  zero  then  it  lies 
outside  the  dataset  defined  by  the  intersection  of  the  associated  grid  line  with  con¬ 
tours,  is  for  a  closed  map  the  point  lies  beyond  the  outermost  contours.  The  value  of 
the  final  grid  at  that  point  should  therefore  have  the  value  ZE  which  is  to  be  assigned 
to  such  regions.  Applying  the  same  procedure  to  open  maps  would  produce  undesirable 
effects,  for  example  all  values  within  the  shaded  area  of  Fig  9b  would  be  the  same. 

The  final  grid  produced  by  combining  those  shown  in  Figs  4  and  5  with  this  method 
is  shown  in  Fig  10. 

4 . 5  Alternative  version 

The  interpolation  algorithm  described  in  this  Report  for  carrying  out  interpolation 
with  contour  data  is  adopted  so  as  to  give  a  'best'  estimate.  However,  if  the  contours 
have  been  digitized  such  that  linear  interpolation  is  appropriate,  an  alternative  version 
of  subroutine  COORDS  employing  linear  interpolation  is  available.  This  subroutine  is 
listed  in  Appendix  D  and  can  be  simply  substituted  for  that  given  in  Appendix  B,  omitting 
subroutines  INTERP  and  CUBIC  which  are  no  longer  required. 

4.6  The  problem  of  interpolation  outside  the  contour  set 

In  the  case  of  closed  maps,  the  user  of  the  program  is  required  to  decide  what 
values  should  be  assigned  to  grid  points  lying  outside  the  outermost  contour.  Two  poss¬ 
ible  options  are: 

(i)  Setting  all  values  to  the  same  value,  say  that  of  the  base  level  or  the 

outermost  contour. 

(ii)  Extrapolating  smoothly  to  some  defined  level. 

The  first  option  is  easy  to  implement  and  is  available  with  this  program,  but  the 
second  is  much  more  difficult.  It  should  be  realised  that  one  cannot  reliably  extrapolate 
beyond  the  outermost  contour  unless  the  nature  of  the  variable  in  the  region  is  known.  It 
may  be  that  the  user  knows  that  the  variable  tends  to  some  base  level  and  so  wants  to 
extrapolate  smoothly  to  this  level  for  'aesthetic'  reasons.  In  this  situation  one  very 
satisfactory  method  is  for  the  user  to  define  a  'false'  contour  beyond  the  outermost  con¬ 
tour,  where  it  is  estimated  that  the  variable  reaches  the  base  level,  define  a  second 
false  contour  (the  shape  of  which  is  not  important)  outside  the  first  and  then  run  the 
program  with  all  values  beyond  the  false  contours  being  set  to  the  base  level.  The 
second  false  contour  ensures  that  along  a  grid  line  the  interpolation  curve  will  have  a 
gradient  of  zero  when  it  reaches  the  base  level  and  that  areas  of  constant  value  do  not 
occur  in  the  region  between  the  first  false  contour  and  the  outermost  contour  of  the  ori¬ 
ginal  data.  (This  latter  point  is  discussed  in  more  detail  in  section  5.)  The  difficulty 
of  trying  to  extrapolate  analytically  is  that  the  contours  defined  by  the  extrapolation 
routines  may  be  far  from  smooth,  which  is  undesirable  if  the  original  map  contours  are 
smooth . 
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If  the  first  of  the  above  options  is  applied  to  open  maps,  then  the  areas  enclosed 
by  a  corner  of  the  map  and  the  contour  closest  to  that  corner,  will  be  areas  of  constant 
height.  Thus,  if  possible  the  user  should  enlarge  the  area  over  which  contours  are 
digitized,  such  that  the  required  area  is  contained  within  the  final  grid  but  does  not 
encompass  any  of  the  areas  of  constant  height  occurring  at  the  corners. 

4 . 7  Checking  the  program 

Two  programs  are  listed  in  Appendix  E  to  enable  a  user  to  test  program  MAPGRID. 

The  program  DRIVER!  will  produce  the  contour  data  that  was  used  in  running  MAPGRID  to 
produce  the  grid  shown  in  Fig  10.  Program  DRIVER2  will  produce  the  same  grid  analytic¬ 
ally,  the  difference  between  the  two  being  that  shown  in  Fig  8. 

5  USING  PROGRAM  MAPGRID 

For  most  contour  maps  that  it  is  anticipated  that  a  user  will  encounter  the  method 
described  in  this  paper  will  work  well.  There  are,  though,  several  situations  where  the 
resultant  grid  is  not  particularly  reliable.  Here  some  examples  of  this  will  be 
discussed  together  with  some  methods  for  overcoming  the  difficulties. 

In  Fig  8  the  greatest  difference  between  the  interpolated  and  analytic  grid  is  5, 
compared  with  a  contour  interval  of  25,  which  generally  would  be  considered  to  be  good 
agreement  since  no  a  'priori  knowledge  of  the  variable  is  assumed,  other  than  it  be 
smooth  and  continuous.  However,  it  is  worth  examining  a  diagonal  cross-section  which 
passes  through  the  centre  of  the  grid.  The  resultant  curve  (which  is  produced  by  pro¬ 
gram  SECTION  described  in  section  6.2)  is  shown  in  Fig  11,  where  it  can  be  seen  that 
there  is  a  slight  'pecularity'  in  the  curve  near  contour  level  200.  The  cross-section 
is  at  45°  to  the  two  sets  of  grid  lines  and  it  is  evident  that  if  the  cross-section  had 
been  taken  parallel  to  a  set  of  grid  lines  the  feature  would  not  occur.  The  reason  for 
the  occurrence  of  this  feature  can  best  be  illustrated  by  reference  to  Fig  12.  For  grid 
line  aa  the  interpolated  values  in  the  range  bb  become  increasingly  unreliable  as 
bb  ,  and  hence  d  ,  increases.  Interpolated  values  are  most  unreliable  for  the  largest 
values  of  d  which  satisfy  the  condition  that  d  be  less  than  r2  ~  r\  >  ^ut  interpo¬ 
lated  values  in  the  interval  bb  are  not  incorporated  in  the  grid  if  the  weighting  of 
the  points  associated  with  the  orthogonal  set  of  grid  lines  is  greater.  The  weighting 
values  associated  with  the  final  grid  shown  in  Fig  10  can  be  seen  in  Fig  13.  The 

weighting  values  indicated  by  1  lie  near  point  X  in  Fig  12.  By  inspection  it  can  be 

seen  that  for  certain  values  of  r^  and  r(  ,  X  does  not  lie  inside  r^  ,  the 
criterion  for  this  being  that  r^/r ^  <  /2  .  This  is  why  no  other  feature  similar  to 
that  which  occurs  at  contour  level  200  exists  on  the  cross-section,  ie 

r 1 75  .  ,  r200  . 

-  =  1.316  whereas  -  ”  1.681 

r200  r225 

where  the  subscripts  indicate  the  contour  heights.  Thus,  when  taking  a  series  of  cuts 

across  a  map  it  is  preferable  to  take  these  parallel  to  a  set  of  grid  lines. 
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Another  difficulty  which  may  be  encountered  is  best  illustrated  by  reference  to 
Fig  l-ta.  In  this  figure  there  are  two  sets  of  closed  contours  indicated  by  A  and  H  , 
set  A  consists  of  three  contours  and  set  B  of  onlv  one  contour.  By  inspection  it 
can  be  seen  that  within  the  contour  of  set  B  the  variable  is  greater  than  I  and  less 
than  the  next  contour  level  of  2,  and  outside  sets  A  and  B  the  variable  is  less 
than  I  but  greater  than  or  equal  to  the  base  level  of  0 .  It  is  not  possible  however  to 
make  any  reliable  estimate  as  to  the  magnitude  of  the  departure  of  the  variable  from  I  it; 
these  cases.  Fig  lAbAc  indicates  the  reliability  of  interpolated  values  obtained  by 
taking  cuts  across  the  map  in  the  x  and  y  directions.  Fig  1 4 cl  shows  the  values 
obtained  after  the  two  grids  have  been  combined,  from  which  it  can  be  seen  that 
reliable  information  as  to  the  shape  of  the  lowest  contour  is  lost.  The  contour  with 
value  2  is  correctly  described  and  all  interpolated  values  within  it  are  considered 
reliable. 

The  lowest  contour  level  and  the  sense,  although  not  the  magnitude,  of  the  inter¬ 
polated  data  can  be  correctly  described  using  the  following  method.  The  contour  map  is 
considered  as  consisting  of  two  separate  contour  maps,  these  being  the  contour  sets 
A  and  B  .  Another  contour  with  a  value  between  0  and  1  is  drawn  outside  the  two  sets 
as  shown  in  Fig  15.  If  program  MAPGRID  is  now  run  (IROUTE  set  to  I,  ZE  set  to  the  base 
level  of  0),  the  resultant  grid  will  be  reliable  for  grid  values  greater  than  or  equal  to 
1.  It  is  suggested  when  drawing  the  false  contours  that  the  closer  the  chosen  contour 
level  of  the  false  contour  is  to  the  value  of  the  outermost  contour  the  closer  the  false 
contour  should  be  drawn  to  that  contour.  In  drawing  a  false  contour  close  to  the  outer¬ 
most  contour  a  user  is  most  likely  to  estimate  accurately  the  behaviour  of  the  variable 
beyond  the  outermost  contour.  This  procedure  of  drawing  false  contours  is  recommended 
even  if  the  contour  map  shown  in  Fig  14a  were  to  consist  only  of  contour  set  A  ,  since 
from  inspection  of  Fig  I4d  certain  parts  of  the  interval  between  contour  levels  I  and  2 
contain  areas  of  constant  height  not  implied  by  the  data. 

The  above  discussion  is  by  no  means  exhaustive,  but  is  merely  included  to  indicate 
some  problems  which  can  be  encountered  and  how  it  may  be  possible  to  overcome  them. 
Clearly  the  user  should  study  any  map  before  using  the  technique  described  in  this 
Report  to  ascertain  whether  it  is  desirable  to  introduce  any  false  contours  or  to  con¬ 
sider  the  map  in  several  parts. 

6  ANCILLARY  PROGRAMS 

Two  further  programs  have  been  written  to  carry  out  operations  on  the  grid  pro¬ 
duced  by  program  MAPGRID.  One  of  the  programs  expands  the  number  of  data  points  in  the 
grid  whilst  the  other  enables  a  cross-section  to  be  taken  through  it. 

6. 1  Program  EXPAND 

One  of  the  requirements  that  the  author  has  is  to  display  the  grid  on  a  TV  moni¬ 
tor  so  as  to  visually  inspect  the  data.  The  number  of  pixels  required  to  fill  the 
screen  has  been  512  *  512,  a  number  considerably  greater  than  that  needed  in  other 
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..'mpulor  pro.  rssi  n^;  ,>t  tin1  grid.  This  program  lias  therefore  been  written  to  expand  (or 
contra.  i  )  t  lu  number  ot  data  points  in  either  or  both  the  x  and  y  directions.  The 
different  expansion  tutors  for  the  two  directions  means  that  the  aspect  ratio  (width/ 
height*  ot  the  grid  i  an  be  altered  so  that  the  resulting  image  appears  correctly  pro¬ 
portioned  on  a  IV  monitor  (normal  aspect  ratio  of  1.33). 

(are  tias  to  he  exercised  when  interpolating  data  that  has  already  been  interpolated. 
In  general  it  is  accepted  that  in  such  instances  it  is  best  to  use  the  same  interpolation 
algorithm  as  was  used  to  produce  the  interpolated  data.  Thus  in  this  program  the  average 
tangents  algorithm  t  section  i.  1.1  was  used,  but  it  omits  the  ’logic'  which  prevents 
data  at  turning  points  reaching  the  next  contour  level  as  information  concerning  the  next 
contour  level  is  not  contained  in  the  grid. 

The  operation  ot  the  program  is  straightforward,  the  user  only  having  to  specify 
the  expansion  factor  repaired  in  each  of  the  directions.  The  number  of  points  down  each 
column  of  the  grid  is  tirst  adjusted  to  the  required  number  and  then  the  procedure  is 
repeated  tor  each  line,  each  line  being  written  to  an  output  file  as  it  is  computed.  A 
desi  r i pt i on  and  a  listing  ot  the  program  is  given  in  Appendix  F. 
b  .  J  Program  Sht.TlON 

A  co mm. n  requirement  is  to  look  at  a  section  across  a  grid.  This  program  enables 
anv  eross-sei  t i on  in  anv  direction  to  be  taken.  Again,  for  consistency,  the  average  tan¬ 
gents  algorithm  is  used  for  interpolation. 

The  method  adopted  is  as  follows.  The  co-ordinates  at  which  interpolated  values 
along  the  cut  are  required  are  calculated  from  the  co-ordinates  of  the  two  ends  of  the 
cut  and  the  number  ot  points  required  in  the  cut.  If  the  number  requested  is  less  than 
two  then  it  is  assume. l  that  points  are  required  at  the  same  interval  as  the  grid  spacing. 
Values  at  the  required  points  (*  in  Fig  16)  are  evaluated  as  follows.  Four  points  (o) 
are  calculated  with  the  average  tangents  algorithm  using  the  four  points  linked  by  the 
dashed  line.  These  four  points,  linked  by  the  dotted  line,  are  then  used  to  obtain  a 
value  at  the  required  point  (*).  This  procedure  is  then  repeated  in  the  orthogonal  direc¬ 
tor!  and  a  second  estimate  obtained.  The  final  value  is  taken  as  the  average  of  the  two 
estimates.  In  practice  it  is  found  that  the  difference  between  the  two  estimates  is 
insignificant  if  the  grid  has  been  produced  by  program  MAPGRID.  With  this  method,  inter¬ 
polated  values  cannot  be  obtained  in  the  end  intervals,  so  values  required  in  these 
intervals  are  initially  ignored  and  only  those  for  which  the  above  method  is  applicable 
are  derived.  Tin  values  required  in  the  end  intervals  are  obtained  from  linear  extrapola¬ 
tion  of  the  data  previously  derived. 

A  description  and  a  listing  of  the  program  is  given  in  Appendix  G. 

7  CONCLUSIONS 

This  Report  has  described  a  technique  for  converting  digitized  contour  data  into 
an  equi spaced  grid  of  points,  and  a  computer  program  MAPGRID  which  performs  this  trans¬ 
formation  has  been  included.  It  should  be  realised  that  the  computer  programs  described 
in  this  Report  are  not  designed  either  for  maximum  efficiency  or  minimum  storage  require¬ 
ments.  They  are  included  so  that  the  reader  can  use  them  to  assess  the  technique  and  are 
written  for  clarity  so  that  they  can  easily  be  modified  if  required.  Some  applications 
of  this  technique  are  to  be  described  in  a  subsequent  report. 


MAPCRID  MAIN  PROGRAM 

A. I  Funct ion 

The  main  program  initialises  constants,  supervises  the  flow,  controls  all  input 
and  output  and  combines  the  data  from  the  two  sets  of  cuts  to  produce  the  final  grid. 

A. 2  Dimensions  of  arrays 

F  l 

.  ^  maximum  number  of  points  in  any  contour 

X  ) 

Y  /  maximum  number  of  points  of  intersection  of  a  set  of  grid  lines  with  the  contours 

Z  I 

A  I 

S  Y  maximum  number  of  contour  crossings  of  any  grid  line 

W  I 

l'  number  of  grid  lines  parallel  to  the  y  axis  (1DIM) 

V  number  of  grid  lines  parallel  to  the  x  axis  (JDIM) 

WT  the  larger  of  IDIM  and  JDIM 

H  number  of  contours  +2 

W \  }  IDIM’  mLM- 

A. 3  Constants  to  be  initialised 

MAX  maximum  number  of  contour  crossings  of  any  grid  line  (»  dimension  of  A,  S  and  W) 

KMAX  maximum  number  of  points  of  intersection  of  a  set  of  grid  lines  with  the  con¬ 

tours  (=  dimension  of  X,  Y  and  Z) 

IDIM  number  of  grid  lines  parallel  to  the  y  axis 

JDIM  number  of  grid  lines  parallel  to  the  x  axis 

IROL'TE  control  parameter:  =  I  for  closed  maps,  *1  for  open  maps 

ZE  value  to  be  assigned  to  points  outside  the  contour  map 

WV,'  grid  spacing 

XST  x  co-ordinate  of  the  origin  of  the  grid 

YST  y  co-ordinate  of  the  origin  of  the  grid 

HMIN  minimum  allowed  value  of  map  data 

HMAX  maximum  allowed  value  of  map  data 

FRACTN  normally  set  at  0.99,  see  section  3.3.2 
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A. 4  Program  units 

CHECK  ensures  that  no  two  adjacent  data  points  are  coincident 
IDENT  reformats  the  contour  data 

COORDS  evaluates  the  co-ordinates  at  which  y  direction  grid  lines  intersect  a 
contour 

NFIX  gives  the  integer  value  of  a  number,  the  value  being  rounded  down  to  the  lower 
integer 

INTERP  calculates  the  co-ordinates  at  which  a  line  parallel  to  the  y  axis  intersects 
a  contour  between  two  points 

CUBIC  finds  the  roots  of  a  third  order  polynomial  in  x  in  a  specified  range 

HEIGHT  sorts  the  contour  values  into  ascending  order 

SELECT  selects  data  points  with  a  given  x  co-ordinate 

CHECK2  ensures  that  no  two  data  points  on  a  grid  line  are  coincident 

ORDER  sorts  values  in  an  array  into  ascending  order 

AVTAN  calculates  equispaced  grid  points  along  a  grid  line 

AVTAN2  performs  average  tangents  interpolation  between  two  data  points 

WEIGHT  generates  weighting  values  associated  with  interpolated  data  points. 

In  addition  the  following  standard  FORTRAN  functions  are  called: 

FLOAT  converts  integer  to  floating  point 

IFIX  integer  value  nearer  to  zero 

ABS  absolute  value 

ATAN  arc  tangent 

COS  cosine 

A. 5  Input/output 

The  contour  data  being  used  by  the  program  is  required  to  be  read  twice.  It  is 
input  to  the  program  via  FORTRAN  channels  5  and  7.  The  output  data  is  written  to  an 
output  file  via  FORTRAN  channel  6. 


\ 
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Appendix  B 

SUBROUTINE  SPECIFICATIONS  FOR  PROGRAM  MAPCRID 
SUBROUTINE  CHECK 

~  This  subroutine  ensures  that  no  two  adjacent  data  points  are 
coincident 

-  SUBROUTINE  CHECK  (F,  C,  IMAX,  NMAX) 

-  NMAX  dimension  for  arrays  F  and  G 

£  |  x  and  y  co-ordinates  of  the  data  points 
IMAX  number  of  data  points 
Subordinate  subprograms  -  None 

Lxp 1 anat ion  -  A  check  is  carried  out  with  each  pair  of  adjacent  data  points 

to  determine  whether  they  have  the  same  co-ordinates.  If  they  do  one  of  the  points  is 
removed  and  IMAX  is  decremented  by  I. 


Summary 

Subroutine  statement 

Input  argument 

Input / output  arguments 


SUBROUTINE  IDENT 


Summary 


-  This  subroutine  modifies  the  format  of  contour  data  so  that 
it  is  in  a  suitable  form  for  subroutine  COORDS 


Subroutine  statement 


Input  arguments 


Input/output  arguments 


-  SUBROUTINE  IDENT  (F,  G,  IMAX,  NMAX,  ITYPE) 

-  NMAX  -  dimension  for  arrays  F  and  G 

ITYPE  -  contour  type:  open  contour  «  +!,  closed  contour  »  0 

-  F  \  x  and  y  co-ordinates  of  data  points  defining  the 
C  f  contour 

IMAX  number  of  data  points  defining  the  contour 


Subordinate  subprograms  -  None 


Exp 1 anat ion 


See  section  A. 2. 
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SUBROUTINE  COORDS 

Summary 

Subroutine  statement 
Input  arguments 


Output  arguments 

Input /output  argument 


vious  and  present  cans  to  tne  subroutine 
since  the  initialization  of  K  in  the  main 
program 

Subordinate  subprograms  -  Function  NFIX 

Subroutine  INTERP 

Explanation  -  Each  pair  of  data  points  between  2  and  IMAX-)  is  taken  in  turn 

and  the  y  direction  grid  lines  which  intersect  the  contour  between  these  two  points  are 
derived.  Subroutine  INTERP  is  called  to  find  the  point  of  intersection  of  each  of  these 
grid  lines  with  the  contour.  For  each  point  of  intersection  K  is  incremented  by  I  and 
the  co-ordinates  are  entered  into  arrays  X,  Y  and  Z.  It  is  possible  that  the  contour  is 
changing  direction  between  the  two  data  points  with  the  result  that  one  or  more  grid 
lines  may  cross  the  contour  in  two  places  (see  Fig  17).  In  this  situation,  subroutine 
INTERP  is  called  to  determine  whether  the  contour  cuts  the  first  grid  line  beyond  the 
point  farthest  from  the  point  at  which  the  contour  changes  direction  (P^  in  Fig  17). 

If  it  does  not,  the  next  pair  of  points  (P^  and  P^  +  ^)  ate  considered.  If  it  does,  then 
the  co-ordinates  of  the  intersection  points  are  entered  into  arrays  X,  Y  and  Z  and  K  is 
suitably  incremented.  (If  the  contour  just  touches  a  grid  line  two  points  of  intersec¬ 
tion  are  still  recorded  although  these  points  are  coincident.)  The  next  farther  grid 
line  is  now  taken  then  the  process  repeated  until  a  grid  line  is  found  that  does  not 
cross  the  contour,  when  the  next  pair  of  points  is  then  considered. 

Only  pairs  of  points  between  2  and  IMAX-I  are  considered  because  of  the  nature  of 
the  interpolation  algorithm  used  in  subroutine  INTERP. 

A  parameter  EPS  is  defined  in  the  subroutine,  and  is  given  as  0.01  *  WF  in  the  pro¬ 
gram  listing.  Thig  parameter  should  be  chosen  so  that  it  is  smaller  than  the  digitizing 
step  of  the  digitizer  used  to  digitize  the  contour  data. 


-  This  subroutine  evaluates  the  co-ordinates  at  which  y  direc¬ 
tion  grid  lines  intersect  a  contour 


-  SUBROUTINE  COORDS  (F,  0,  H,  X,  Y,  Z,  IMAX,  KMAX,  WF,  K) 
r  ^  arrays  containing  the  x  and  v  co-ordinates  of 


-  F 

G 


arrays  containing  the  x  and 
contour  data  points 


H 


contour  height 


IMAX  number  of  data  points  defining  the  contour 


KMAX 

WF 


-  X 
Y 
Z 

-  K 


dimensions  of  arrays  X,  Y  and  Z 
grid  spacing 

arrays  containing  the  x  and  y  co-ordinates  and  the 
height  of  intersection  points  between  the  grid  lines 
and  the  contour 

On  input  :  number  of  intersection  points  found  by  pre¬ 
vious  calls  to  the  subroutine  since  the 
initialization  of  K  in  the  main  program 

On  output  :  number  of  intersection  points  found  by  pre- 
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FUNCTION  NFIX 

Summary  -  This  function  gives  the  integer  value  of  a  number,  the  value 

being  rounded  down  to  the  lower  integer 

Function  statement  -  FUNCTION  NFIX(FF) 

Input  argument  -  FF  floating  point  number 

Output  function  -  NFIX  required  integer  value 

Subordinate  subprograms  -  FORTRAN  function  IFIX 

Explanation  -  A  call  is  made  to  a  computer  library  function  IFIX  which  gives 

the  integer  value  of  a  number,  the  value  being  the  nearer  integer  to  zero.  If  the 
input  value  parameter  has  a  value  which  is  not  an  exact  integer,  the  result  of  IFIX  has 
1  subtracted  from  it  if  the  input  parameter  to  IFIX  is  negative.  The  user  is  advised 
that  his  compiler  may  have  a  library  function  equivalent  to  NFIX. 


SUBROUTINE  INTERP 


Summary 


Subroutine  statement 


Input  arguments 


Output  arguments 


-  This  subroutine  calculates  the  co-ordinates  at  which  a  line 
X  =  XX  intersects  a  contour  between  points  (X2,  Y2)  and 
(X3,  Y3) 

-  SUBROUTINE  INTERP  (XI,  Y1 ,  X2,  Y2,  X3,  Y3,  X4,  Y4,  XX,  YY, 

YYS,  ITEST) 

-XI  V 
YI  1 
X2  1 

Y2  lx  and  y  co-ordinates  of  the  points  used  by  the 
X3  I  interpolation  algorithm 
Y3  | 

X4  I 
Y4  / 

XX  x  value  at  which  an  interpolated  y  value  is  required 

-  YY  I 

YYg  (  interpolated  y  values 

ITEST  number  of  interpolated  y  values 


Subordinate  subprograms  -  SUBROUTINE  CUBIC 

Explanat ion  -  The  interpolated  values  of  y  are  found  using  the  average 

tangents  algorithm  described  in  section  3.3.3.  The  equation 

3  2 

x  =  as  +  bs  +  cs+d 
X  X  X  x 


is  solved  by  a  call  to  subroutine  CUBIC. 


■v 
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SUBROUTINE  CUBIC 


Summary 


Subroutine  statement 


Input  arguments 


Output  arguments 


-  This  subroutine  finds  the  roots  of  a  third  order  polynomial  in 
x  in  the  range  0  <  x  <  x 

max 

-  SUBROUTINE  CUBIC  (A,  B,  C,  D,  XI,  X2,  X3,  XMAX,  ITEST) 

3 

-  A  coefficient  of  the  x  term 

2 

B  coefficient  of  the  x  term 

C  coefficient  of  the  x  term 

D  constant  term 

XMAX  maximum  value  of  x  for  which  roots  are  required 

-  XI  first  root 

X2  second  root 

X3  third  root 


ITEST  number  of  roots 

Subordinate  subprograms  -  FORTRAN  functions  ABS 

ATAN 

COS 


-  This  subroutine  gives  the  real  roots  of  the  equation 
ax3  +  bx2  +  cx  +  d  =  0  .  A  check  is  first  carried  out  to  determine  whether  the 


Explanation 
f(x) 

contribution  of  the  x  term  is  significant;  the  criteria  adopted  being  that 

lOOax  >  b  .  If  this  criterion  is  not  satisfied,  the  equation  is  treated  as  a  quadra- 
max 

tic.  For  this  case  a  check  is  carried  out  to  determine  whether  the  contribution  of  the 
x2  term  is  significant;  the  criterion  adopted  being  that  'OObx  >  c  .  If  this 
criterion  is  not  satisfied,  the  polynomial  is  reduced  to  a  linear  equation.  If 


lOOax  >  b  but 
max 


I04ax2 


<  c  the  polynomial  is  again  reduced  to  a  linear  equation. 


The  roots  of  the  equation  for  these  three  cases  are  as  follows: 


0) 


3  2 

Cubic  equation  f (x)  =  ax  +  bx  +  cx  +  d  =  0 


Q  « 


3c  -  b  /a 
9a 


R  = 


9bc/a  -  27d  -  2b3/a2 
54  a 


3  2 

If  Q  +  R  >  0  ,  then 


■ c°*D  "■"'(zb)]  ■  & 


3  2 

If  Q  +  R  >0 


3/r  +  /o3  +  R2 


T  = 


'h- 


/?  ♦  r2 


and 


-  S  +  T  -  — 
3a 


i 

i 

I 


110 


'V 
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Having  found  the  first  root  the  other  two  can  now  be  found. 


(2)  Quadratic  equation  f(x) 


if  -  Tp  <  0 

if  ~  ^  0 


2 

The  roots  are  real  only  if  c  4bd  4-  0  . 
(3)  Linear  equation  f (x)  =  cx  +  d  =  0 


Each  of  the  roots  is  tested  to  determine  whether  it  lies  in  the  range 

0  <  x  <  x  .  ITEST  indicates  the  number  of  roots  that  satisfy  this  criteria.  If 
max 

there  is  only  one  root  it  is  returned  by  XI  and  if  there  are  two  they  are  returned  by 
XI  and  X2  . 
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SUBROUTINE  HEIGHT 

Summary  -  This  subroutine  sorts  the  contour  values  into  ascending  order. 

The  first  and  last  values  are  contour  levels  within  which 
interpolated  values  are  to  be  restricted. 

Subroutine  statement  -  SUBROUTINE  HEIGHT  (H,  M,  IH,  HMIN,  UMAX) 


Input  arguments 


Input/output  arguments 


-  M  number  of  contours  used  in  program  +2 

HMIN  minimum  allowed  value  of  map  data 
HMAX  maximum  allowed  value  of  map  data 

-  H  On  input  :  array  containing  all  the  height  values  of 

all  the  contours  used  in  the  program 


On  output  :  array  containing  the  values  of  all  the 

different  contour  heights  plus  HMIN  and  HMAX 


IH  On  input  :  number  of  contours  used  in  the  program 
On  output  :  number  of  different  contour  heights  +2 


Subordinate  subprograms  -  None 


Explanation  -  HMIN  and  HMAX  are  entered  into  array  H  and  then  the  values 

of  H  are  sorted  into  ascending  order  using  a  bubble  sort  (see  explanation  under  SUB¬ 
ROUTINE  ORDER).  A  test  is  then  carried  out  with  each  pair  of  values  of  H  to  see  if 
any  are  equal.  If  the  two  values  are  found  to  be  equal,  one  of  the  values  is  removed 
from  the  array  and  IH  ,  denoting  the  values  of  interest  in  the  array,  is  decremented 
by  1 . 


Summary 


SUBROUTINE  SELECT 

-  This  subroutine  selects  data  points  with  a  given  x 
co-ordinate 


Subroutine  statement 
Input  arguments 


Output  arguments 


-  SUBROUTINE  SELECT  (X,  Y,  Z,  KMAX,  U,  V,  W,  MAX,  J,  WW) 


-  X 
Y 
Z 


KMAX 

U 

MAX 

WW 

V 

W 

J 


arrays  containing  the  x,  y  and  z  co-ordinates  of  the 
intersection  points  between  grid  lines  and  contours 

number  of  points  of  intersection 

value  of  x  for  which  data  points  are  required 

dimensions  for  arrays  V  and  1' 

grid  spacing 

array  containing  v  values  with  the  x  value  U 

array  containing  z  values  with  the  x  value  U 

number  of  data  points  with  x  value  U 


Subordinate  subprograms  -  None 

Explanation  -  Each  point  of  intersection  is  tested  to  determine  whether  the 

x  value  is  equal  to  U  .  If  it  is,  J  is  incremented  by  I  and  the  corresponding  y 
and  z  values  are  entered  into  arrays  V  and  W  respectively. 
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SUBROUTINE  ORDER 

Summary  -  This  subroutine  orders  values  in  an  array  into  ascending  order 

Subroutine  statement  -  SUBROUTINE  ORDER  (V,  W,  LMAX) 

Input  argument  -  LMAX  number  of  values  to  be  sorted 

Input/output  arguments  -  V  array  of  values  to  be  sorted 

W  array  of  values  associated  with  W 

Subordinate  subprograms  -  None 

Expl anat ion  -  The  method  used  is  the  bubble  sort.  This  involves  taking  each 

pair  of  values  of  V  in  turn  and  checking  to  determine  whether  or  not  they  are  in  the 
correct  order.  If  they  are  not  they  are  interchanged  (together  with  the  associated 
values  of  W  )  and  a  counter  J  is  incremented.  If  J  is  non-zero  (indicating  that  at 
least  one  pair  of  values  has  been  interchanged)  the  process  is  repeated.  When  J  is 
zero  all  the  values  are  in  the  correct  order. 

SUBROUTINE  CHECK2 

Summary  -  This  subroutine  ensures  that  in  the  situation  where  a  grid  line 

just  touches  a  contour,  the  two  data  points  representing  the 
points  of  intersection  of  the  grid  line  with  the  contour  are 
not  coincident. 

Subroutine  statement  -  SUBROUTINE  CHECK2  (V,  LMAX,  WW) 

Input  arguments  -  LMAX  number  of  data  points 

WW  grid  spacing 

Input/output  arguments  -  Y  array  containing  y  values  of  data  points 

Subordinai-e  subprograms  -  FORTRAN  function  ABS 

Explanation  -  If  a  grid  line  just  touches  a  contour  subroutine  COORDS  outputs 

two  points  with  the  same  co-ordinates,  rather  than  one  point,  which  will  cause  subsequent 
interpolation  algorithms  to  fail  unless  corrected.  This  subroutine  checks  each  adjacent 
pair  of  points  to  determine  whether  they  have  the  same  y  values.  If  they  do  the  y 
values  of  these  points,  y^  ^  and  y.  ,  are  then  replaced  by  y.  ^  -  E(  and  y.  +  E^ 
where 

Lj  =  |yi_,  ~  y  £_2  i /5  unless  E  ^  >  E  or  y^_^  c*oes  not  exist,  when  E(  =  E 

e2  -  lyi+i  -  yj/s  unless  E^  >  E  or  y.+^  does  not  exist,  when  E^  *  E 


and  E 


0.01  x  WW  . 


Appendix  B 


18 


SUBROUTINE  AVTAN 

Summary  -  This  subroutine  calculates  equispaced  data  points  along  a  grid 

line 

Subroutine  statement  -  SUBROUTINE  AVTAN  (A,  Y,  Z,  N,  YY,  ZZ,  JDIM,  H,  IH,  FRACTN , 

I ROUTE,  ZE) 


Input  arguments 

-  Y 

array  containing  y  values  of  data  points 

Z 

array  containing  z  values  of  data  points 

N 

number  of  data  points 

YY  ^ 

array  containing  y  values  at  which  interpolated 
values  of  z  are  required 

JDIM 

number  of  interpolated  values  required 

H  1 

IH  | 

FRACTN 

|  used  in  SUBROUTINE  AVTAN2  (see  subroutine  AVTAN2) 

IROUTE 

routing  parameter  (see  explanation) 

ZE 

value  of  z  for  points  outside  the  contour  map 

Output  argument 

-  ZZ 

(  array  containing  interpolated  values  corresponding 
{  to  YY 

-  A 

array  containing  the  average  tangent  at  each  of  the 
data  points 

Subordinate  subprograms  -  SUBROUTINE  AVTAN2 

Explanation  -  The  average  tangents  at  each  of  the  data  points  (y^,  z^)  are 

first  calculated.  For  each  value  of  y  at  which  interpolated  values  of  z  are 
required  (yy)  the  data  points  are  found  which  satisfy  the  criterion  y^  <  yy  <  y^+j  ■ 

A  call  is  then  made  to  subroutine  AVTAN 2  to  obtain  the  interpolated  value  of  z  at  yy  . 
If  the  value  of  yy  lies  outside  the  range  of  the  data  points  then  one  of  two  alterna¬ 
tives  is  taken  depending  upon  the  value  of  the  routing  parameter  IROUTE. 

IROUTE  *  1  The  interpolated  value  is  set  to  the  value  of  z  corresponding  to 
the  y  value  closest  to  yy  . 

IROUTE  =  I  The  interpolated  value  is  set  to  ZE  ,  a  value  specified  by  the  user 
in  the  main  program. 

It  should  be  noted  that  in  arrays  YY  and  Y  the  values  must  be  arranged  in 
ascending  order  and  no  two  values  of  Y  must  be  the  same.  The  latter  condition  is 
ensured  by  a  call  to  subroutine  CHECK2  prior  to  the  call  to  this  subroutine. 
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Summarv 


SUBROUTINE  AVTAN2 

-  This  subroutine  performs  average  tangents  interpolation 
between  two  data  points 


Subroutine  statement 


Input  arguments 


Output  argument 
Subordinate  subprograms 


-  SUBROUTINE  AVTAN2  (Yl,  Y2,  Zl,  Z2,  A1 ,  A2,  YY ,  ZZ,  H,  IH, 


FRACTN) 

-Y1  j 
Zl 

A!  J 

[  y  and  z  values  and  average  tangent 
|  data  point 

at  the  first 

Y2  j 
Z2 

A2  ! 

1  y  and  z  values  and  average  tangent 

(  data  point 

at  the  second 

YY  j 

1  y  value  at  which  an  interpolated  z 

[  required 

value  is 

H  | 

(  array  containing  the  values  of  all  the  different 
contour  heights  plus  minimum  and  maximum  allowed 
(  values  of  map  data 

IH 

number  of  values  in  H 

FRACTN  set  in  main  program,  normally  to  0.99 

-  ZZ 

interpolated  z  value 

-  None 

Explanation  -  The  interpolated  value  of  z  is  found  using  the  average 

tangents  algorithm  described  in  section  3.3.1.  The  algorithm  has  been  modified  (see 
section  3.3.2)  so  that  at  turning  points  the  interpolated  value  is  not  allowed  to 
reach  the  next  contour  level  and  in  other  circumstances  is  constrained  to  remain 
between  Zl  and  Z2  . 
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SUBROUTINE  WEIGHT 

Summary  -  This  subroutine  generates  weighting  values  associated  with  the 

interpolated  data  points  along  a  grid  line 

Subroutine  statement  -  SUBROUTINE  WEIGHT  (Y,  YY,  U,  N,  JDIM,  Z,  A) 


Input  arguments 


Output  argument 


-  Y  array  containing  y  values  of  data  points 

Y'Y  array  containing  y  values  of  interpolated  data  points 

N  number  of  data  points 

JDIM  number  of  interpolated  data  points 

Z  array  containing  z  values  of  data  points 

i  array  containing  the  average  tangents  at  each  of  the 
\data  points 

-  W  array  containing  weighting  values 


Subordinate  subprograms  -  None 


Explanation  -  for  each  interpolated  point  (yy)  data  points  are  found  which 

satisfy  the  criterion  y  <  yy  <  y^+1  .  The  weighting  value  associated  with  yy  is 
given  by: 


Iyi+1  -  V 


tf  y,  =  yy 


If  A.  =  Ai+|  =  0  and  Z.  =  Z.+, 


W  =  -1 


If  yy  lies  outside  the  range  of  y  W  =  0  . 


It  should  be  noted  that  the  weighting  value  is  not  important  when  yy  -  y^  or  y-  +  |  • 


I 
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C  ****41  PROGRAM  MAPGRID  ***** 


PAGE  8001 


c  *****  PROGRAM  HAPCRID  ***** 

C  _  . 

C  ****•***••*****•***••*••***•♦•**•**•***•***•••***•••••••••*•••••••••• 

C  THIS  PROGRAM  COMVERTS  CONVERTS  A  CONTOUR  MAP  INTO  A  UNIFORMLY  SPACED  GRID 

C  OF  POINTS 

C 

C  CONTOUR  FITTING  ALGORITHM i  AVERAGE  TANGENTS 

C 

C  ,,*,*****•**••.•••*••*•••****••***••*••••*••••*•••*••••••*•••••••*•• •••*• 

DIMENSION  F(5B)/C<30)jX<10BB)«Y<  1BBB  >,  2(  1BBB>,  A<  38 ) / 8<  SB  ) 

DIMENSION  U<39>,  V<  59  ),T<39  i.  WT<  39  >.  IK  SB),  M(  SB  > 

DIMENSION  A  A(  39  /  59  >.  W  A<  39  <  39  ) 

INITIALIZATION  OF  PARAMETERS 

MAX-30 
KMAX-1BBB 
IDIM-59 
JDIM-39 
I  ROUTE-! 

ZE-0  0 
ww»i . a 
XST  — 29.B 
YST  — 29  .  B 
HMIN-8  B 
HMAX-23B .8 
FRACTH-0  99 
DO  IB  I-l.IDIM 
IB  U< I  )«FLOAT (  1-1 )*UU 
DO  20  J-l.JDIN 
2B  VC J  )-FLOAT<  J-l )*UU 

GENERATE  CRID  IN  F  DIRECTION 


M-B 

KREC-B 
218  M*M ♦ 1 

READ  (  3<  51  <  END-23B >  H(  N ). I  MAX. I  TYPE 

DO  22B  1-1. IMAX 

READ  ( 5. 32  )  F< I  )*G< I  ) 

FC I )«F< I >-XST 
228  CC I )*G< I  )- YST 
NMAX-INAX 

CALL  CHECK  (FiGj IN  AX  (  NMAX  ) 

NMAX-INAX+2 

CALL  IDENT  < F , C . IMAX , NMAX ,  IT YPE  ) 

CALL  C00RD8  (  F, G.HCM >,X» Y. 2. IMAX, KMAX# 88. KREC  ) 

GO  TO  218 
238  H-M*l 

CALL  HEIGHT  <H<M,IH»HHIN»HMAX> 

DO  248  I-t  <  ID  IN 

CALL  SELECT  < X. Y . Z , KREC , U< I ) . S , V. MAX , LMAX , MW > 

IF  (LMAX  LT.l)  GO  TO  2«B 
CALL  ORDER  <8, W, LMAX) 

CALL  CHECK2  <8,LMAX,NU> 

CALL  AVTAN  < A , S , W, LM AX . V, T , JDI M , M , I H # FRACTN , I  ROUTE , ZE ) 
CALL  WEIGHT  < S, V,WT, LMAX. JOIN, W, A > 

GO  TO  2SB 


non  o  o  o  non 


V 
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26B  DO  271  J-l.JDIH 
HT< J  >--2  I 
27fl  T<  J  )-2E 
291  00  241  J-l.JDIN 
AA(  I.  J  >-T<  J  ) 

241  WA(  I.  J  )-WT(  J> 

GENERATE  GRID  IN  C  DIRECTION 

KREC-B 

311  READ  (7.31. END«33B >  HH. IHAX. ITVPE 
DO  321  I-l . IRAN 
READ  (7.32)  FCI  ).C(I  ) 

F(I >-F( I >-KBT 
321  G(I  )-G(I >-Y8T 

II  IIAWb  t  maw 

CALL  CHECK  ( F . G . IHAX < HR AN  ) 

NHAX- I NAX+2 

CALL  IDENT  (G.F. IRAN, NHAX. ITYPE > 

CALL  COORD8  < C. F . HH. X . Y . 2 . INAX . KHAX , UW . KREC > 

CO  TO  311 

331  WRITE  (6.33)  IDIR.JDIH 
DO  341  J-l.JDIH 

CALL  SELECT  ( X, Y . 2 . KREC . V( J ) . S . U. RAX , L RAX . HU > 

IF  (LHAX  LT .  1  >  GO  TO  361 
CALL  ORDER  (S.U.LHAX) 

CALL  CHECK2  (S.LNAX.UU) 

CALL  AVTAN  ( A. S . U. LHAX. U. T . IDI R , H . I K. FRACTN . I  ROUTE . 2E  ) 
CALL  HEIGHT  ( 8. U . HT. LHAX. ID1H. H . A > 

GO  TO  331 

361  DO  371  I-l .  ID  IR 
UT( I  )*-2 . ■ 

371  T ( I  )-2E 

33B  DO  341  I-l .  IDIR 

IF  (HAd  .  J  )  EG  .  WT(  I  )  >  AA< I.J  )-<AA(I.J)*T( I  >  >/2  .  B 
IF  (HAd.J  ).EQ.WT(  I))  GO  TO  34B 
IF  ( IROUTE  HE  .  1  )  GO  TO  343 
IF  (WA< I ,J  )  EG. -2. E>  GO  TO  34B 
IF  (HT(  I  ).  EQ.  -2  .  ■>  AA(I.J)-Td) 

IF  (HTd). EG. -2  B>  CO  TO  341 
IF  (HAd.J  ).EQ  B  B)  GO  TO  348 
IF  (HTd  ).  EO.  B.  fl  )  AA(  I .  J  )-T<  I  ) 

IF  (HTd  )  EQ.fl  B)  GO  TO  34B 
IF  (HT(I).GT.  UA(  I.J))  AA(I.J)-Td) 

GO  TO  34B 

343  IF  (UA(I.J). EG. -1.8. AND. HT(I  )  LE  B  B>  GO  TO  34B 
IF  (HTd  )  EG  -1  .  B.  AND  HA(  I.J  >  LE.  B  B  >  GO  TO  346 
IF  (HTd  ).GT.UA(  I,  J))  AA(I.J)-Td) 

GO  TO  34B 
346  AA(  I.  J  )- T(  I  ) 

341  CONTINUE 

OUTPUT  FINAL  GRID 

DO  4BB  J-l.JDIH 

4BB  WRITE  (6.34)  (  AAd  .  JD  IN*  1  -  J  ) .  I  ■  1 .  ID  I H  > 

INPUT/OUTPUT  FORMATS 
31  FORMAT  (F6  1. IS.  12) 


noon  o  on  o  o  on  non 
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c  MM*  PROGRAM  MAPGRID 


PACE  0003 


92  FORMAT  < 2F8  3 > 

93  FORMAT  <  13. IX. 13  ) 

94  FORMAT  < 12FC  1 ) 

STOP 

END 

SUBROUTINE  CHECK  < F, G, IMAX. NMAX  ) 

THIS  SUBROUTIHE  ENSURES  THAT  NO  TMO  ADJACENT  DATA  POINTS  ARE  COINCIDENT 

DIMENSION  F<  NNAX  )j  C<  NMAX  > 

1*1 

IB  I -I ♦ 1 

IF  <  I  GT  IMAM  >  RETURN 
IF  <F<  I).NE  F( 1-1 )>  CO  TO  lfl 
IF  (G< I)  ME  C<I-1)>  CO  TO  IB 
DO  28  J*I, IMAX 
F< J-l >-F< J  ) 

2B  C< J-l >-G< J ) 

I  MAX* I  MAX- 1 
I -I -l 
GO  TO  IB 
END 

SUBROUTINE  IOEMT  ( F , G , I  MAX . NMAX  < I  TYPE ) 

THIS  SUBROUTINE  NOOIFIES  THE  INPUT  DATA  FOR  SUBROUTINE  COORDS  DEPENDING 
UPON  THE  CONTOUR  IDENTIFICATION  GIVEN  BY  ITYPE 
ITYPE--U  OPEN  CONTOUR 
ITYPE-  B  CLOSED  CONTOUR 
I TYPE--1  UNDEFINED 

DIMENSION  F<NNAX  ), G< NMAX) 

IF  < I T YPE  >  10,20/ 3B 
IB  RETURN 
2B  F< IMAX+1  )-F< 2  ) 

G< I  MAX* 1  )»G( 2  ) 

F  < I  MAX  *2  )-F<  3  > 

GC IMAX+2 >-G( 3  > 

I  MAX* I MAX*2 
RETURN 

3B  F  C I  MAX ♦ 1  >»2  0*F< IMAX  )-F< IMAX-1 > 

G( I  MAX* 1 >-2.B*G< IMAX  )-G< IMAX-1  ) 

IMAX-IMAX+2 

DO  4B  1-2. INAX 

F< IMAX+2-I )-F< IMAX  +  1-1  ) 

4B  GCIMAX+2-I  )*CCIMAX*1-I> 

F < l >-2B*F< 2>-F<  3) 

G< 1  )-2  B*G( 2)-G<  3  ) 

RETURN 

END 

SUBROUTINE  COORDS  ( F , C, H, X , Y . Z , IMAX . KMAX , MF . K  ) 

THIS  SUBROUTINE  EVALUATES  THE  CO-ORDINATES  AT  UHICH  THE  GRID  LINES 
INTERSECT  A  CONTOUR.  AVERAGE  TANGENTS  INTERPOLATION  IS  EMPLOYED. 

DIMENSION  F< IMAX ), G< IHAX>.X< KM AX > , Y< KMAX  > , 2( KM AX  ) 

EP8-ABS<  B. 01*MF  > 

1-3 

I  FI  -NF  1X(  F  (  1-1  >/'BF  ) 

IF  (A88(  FLOAT! IF1 )*UF-F< 1-1 )  ). GT. EPS  >  GO  TO  9 
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PACE  BBB4 


X-X  +  l 
X(X  )-F< I  ) 

Y(K  )*G<  I  ) 

Z(X  )-H 
3  I -I  - 1 
BB  I -I ♦ 1 

IF  < I  CT  I  MAX  - 1  )  RETURN 
IF1»NF IX<  F< 1-1 >/MF  ) 

I  F2*NF  IX<  F(  D/WF  > 

FF»(F<  IM  )-F<  !))•<  FC  1-1  >-F<  1-2  )  ) 

IF  (FF  LT.I.l)  CO  TO  SB 
NF-IF2-IF1 
IF  <NF)  4B / SB / 6B 
6B  DO  2B  M-l.NF 
X  *K  ♦  1 

X(K  )-FLOAT( 1F1 >*BF*FLOAT<  N  >*MF 

IF  ( ABS< X< X  )-F< I  ))  GT  EPS  )  GO  TO  IB 

Y<K»C<I  ) 

GO  TO  2fl 

IB  CALL  INTERP(F(I-2)»G(  1-2  )  <  F<  1-1  >,  G(  I  -1  >»  F(  I  >»G(  I  >«  FC  I  ♦!  ) 

•  >  G<  !♦  1  ),  X(  X  >.  Y(  X  )  /  YY6 » I  TEST  ) 

2B  Z(X)-H 

CO  TO  SB 
4B  NF--HF 

IF  (AB8<FL0AT(IF2)*«F-F(I  )>  LT  EPS)  NF-NF+1 
DO  3B  M-l.NF 

XX"FLOAT< IF  1 >«BF-FLOAT<  N- 1 >*MF 
IF  CABS' XX-F< 1-1  ))  LT  EPS  )  GO  TO  3B 
X»K  ♦! 

X(X  >-XX 

IF  <ABS(X<  X  )-F<  I  ))  GT  EPS  )  CO  TO  78 
Y<K  )»C<  I  > 

Z(K  )»H 
GO  TO  3B 

7B  CALL  INTERPCFC  I-2>,GU-2),F<  1-1  >.C<  I-t  >.F(  l  >»  G(  I  >,  F<  IM  > 

•  G(  !»1  )/  X<  X  )/ Y<  X  >,  YYS,  I  TEST) 

Z<X  )-H 

3B  CONTINUE 
CO  TO  88 

5B  XX"FLOAT( IF1  )*NF 

IF  (ASS' XX-F< I >  )  GT. EPS  )  CO  TO  S3 
X  «X  ♦  1 
X(X  )-F<  I  ) 

Y(K  >-G<  I  ) 

ZCX  )-H 

35  IF  (FF  CE  B  0 >  CO  TO  BB 

N-B 

ISIGN-t 

IF  << F( !♦! )+F< 1-2)  )-< F< I >*F< I- 1  )>  )  1BB.8B.9B 
SB  ISIGN--1 
N»1 

IF  < I F2 . GT . IF  1 )  1F1-IF2 
GO  TO  1  IB 

IBB  IF  <  IF2.LT  .  IFI  )  IF1MF2 
118  H-N+1*ISIGN 

XX*FLOAT<  IF1+N)*UF 

CALL  INTERP(F<  I-2)/G(  I-2),F<I-1).C<I-1>.F<I>,G(I>.F<I*1) 

•  .  G<  !♦! )»  XX, YY1, YY2/ITEST) 

IF  (ITE8T  EQ  B)  GO  TO  SB 
X-X»l 


uuuu  u  u  u  u  u  u  u  u  o  u  u 
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c  •••••  PROGRAM  MAPGRIO 
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X(  K >-XX 

Y  (  K  >■  Y  Y 1 
2(  K  >-H 

IF  (ITEST  EO  1>  GO  TO  1  IB 
K  «K  ♦  1 
X<K  >-XX 

Y  <  X  T-YY2 
2(K  >»M 

CO  TO  1 1  ■ 

END 

FUNCTION  NF IX<  FF  ) 

THIS  FUNCTION  GIVES  THE  INTEGER  VALUE  OF  A  NUMBER.  THE  VALUE  BEING 
ROUNDED  DOWN  TO  THE  LONER  INTEGER 

NFIX-IFIX(FF) 

IF  (FLOAT(NFIX)  ES  FF)  RETURN 
IF  (FF  LT  ■  ■>  NFIX-NFIX-1 
RETURN 
END 

SUBROUTINE  INTER P  <  XI .Yl, X2. Y2 , X3.Y3.X4. Y4 , XX. YY.YYS.  ITEST  ) 

THIS  SUBROUTINE  CALCULATES  THE  CO-ORDINATES  AT  MHICH  A  LIME  X-XX 
INTERSECTS  A  CONTOUR  BETWEEN  POINTS  (X2,Y2>  AND  (X3.Y3) 

Dl-(( X2-X1  )**2»( Y2-Y1 >**2)*»B  5 
D2-(  <  X3-X2 )**2+(  Y3-Y2  )**2 )**B  5 
D3-((  X4-X3  )**2»(  Y4-Y3  )*«*2  )**B  3 
AXl-(  (  X2-X1  >/Dl*(X3-X2>/»D2  )/2  B 
AX2-((X3-X2  >/02*(X4-X3)/03  )/2  .  B 
A  Y  l  -  (  (  Y2-Y1  )/'Dl*(Y3-Y2)/'D2  )/2  B 
AY2-( (  Y3-Y2  VD2*<  Y4- Y3  )/D3  )/2.  B 
DX-X3-X2 
DY-Y3-Y2 

AX»(AXl*D2+AX2*D2-2  B*DX)/D2**3 
BX-( 3  B*DX-AX2*D2-2  B*AX1*D2  >/0  2**2 
CX-AXt 

AY“(AYl*D2*AY2*D2-2  B*DY  )/02**3 
8Y-( 3 . B*DY- AY 2*02-2  B*A Y 1 *02  )/D 2* *2 
CY-AY1 
DX-X2-XX 

CALL  CUBIC  (AX, BX.  CX , DX , D S 1 , DS 2 , D S3 . D2  , ITEST) 

IF  (ITEST  EO  B>  RETURN 

YY-AY«DS1**3»BY*DS1**2*CY*DS 1»Y2 

IF  (ITEST  EO. 1)  RETURN 

YYS»AY*D82**3*BY*0S2**2*CY*DS2+Y2 

I  TEST-2 

RETURN 

END 

SUBROUTINE  CUBIC  ( A, 8 , C , D , X I , X 2 , X 3 , X MAX ,  I TE ST  ) 

THIS  SUBROUTINE  FINDS  THE  ROOTS  OF  A  THIRD  ORDER  POLYNOMIAL  IN  THE 
RANGE  B  B  LT  X  IT  XHAX 

IF  (ABS( A-XMAX-1BB  B  )  LT  ABS(B>>  GO  TO  3B 

IF  ( ABS( A*XMAX*«2«lBflflB  B)  LT  ABS(C>>  GO  TO  3B 

SOLVE  THIRD  ORDER  POLYNOMIAL 

Q-(  3  B*C-(  B-*2  )/P  ')/(.  9  B*A  ) 
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R-< 9  8*B*C/A-2? . 0*0-2  B-B--3/A  —  2  )/( 54  .8* A) 

IF  (8**3*R**2  LT  .8.8)  GO  TO  5 
S*R+<  0**3*R**2)**8 .5 
SIGH-1  8 

IF  <  S . LT . 8 . 8 )  SIGN  — 1  .8 
S-SICN*(ABS( $))•*< 1.8/38) 

T-R-<  Q**3+R**2>**0 .5 

SIGN- 1  .8 

IF  (T.LT.B.B)  SIGH— 1  .8 
T-SIGH*(ABS<T  )>«•< 1. 8/3.8) 

X 1-S+  T-B/< 3  . B*A  ) 

GO  TO  18 

3  DUM-ATAN( < ABS( -0**3/R**2-l  B))**B  5) 

IF  (R  LT.B.8)  DUN-3. 141592654-DUM 
X 1-2  B*( -8  )**8 . 3*C0SC DUN/ 3  8  )-B/<  3  B-A  ) 

IB  I  TEST- 1 
BB-X1 *B/A 

CC-X1 *< X 1  +  B/A  )*C/A 
TEST-I BB/2  8>**2-CC 
IF  (TEST  LT  .8.8)  GO  TO  25 
IF  ( -BB/2  B.LT.B  8)  GO  TO  15 
X3--BB/2  B+TEST--B .5 
GO  TO  28 

13  X3—BB/2  fl-TEST  **8 . 3 
28  X2-CC/X3 
I  TEST -3 
23  GO  TO  43 

SOLVE  SECOND  ORDER  POLYNOMIAL 

38  IF  ( ABS( B-XHAX-118  8)  LT.ABS(C))  CO  TO  48 
TEST-C  —  2-4  B-B-D 
IF  (TEST  CE  8  I)  CO  TO  33 
I  TEST -8 
RETURN 

33  Xl-(-C-TEST**B  5  )/(2 . B*B) 

X2-( -C  +  TEST--B  3  )/< 2.8*8) 

I  TEST -2 
GO  TO  45 

SOLVE  FIRST  ORDER  POLYNOMIAL 

48  XI— D/C 
I  TEST  - 1 

WHICH  ROOTS  LIE  IN  THE  RANGE  8  8  LT  X  LT  XHAX 

43  IF  (XI  GT  .B  8  AND  XI  LT  XHAX  >  CO  TO  65 
I TE  ST- I  TEST-1 
IF  < I  TEST-1  )  38.33/68 
38  RETURN 
33  X1-X2 
CO  TO  43 
68  X1-X2 
X2-X3 
GO  TO  43 

63  IF  (  ITEST  EM)  RETURN 

78  IF  (X2  GT  8  8  AND  X2  LT  XMAX >  CO  TO  73 
I TEST-1TEST-1 
IF  ( I  TEST  EB  1  )  RETURN 
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X2-X3 
GO  TO  7B 

79  IF  ( I  TEST  EB. 2)  RETURN 

IF  <X3  GT.B  B. AN0.X3.LT. XMAX  >  RETURN 

I  TEST* 1  TEST-1 

RETURN 

END 

SUBROUTINE  HE1CNT  < H , M, IN . HN IN , HNAX > 

C 

C  THIS  SUBROUTINE  SORTS  THE  CONTOUR  VALUES  INTO  ASCENDl NC  ORDER.  THE  FIR8T 

C  AND  LAST  VALUES  ARE  CONTOUR  LEVEL8  WITHIN  WHICH  INTERPOLATED  VALUES  ARE 

C  TO  BE  RESTRICTED 
C 

DIMENSION  H<  H ) 

H( H-l  )-HHIH 
H< M  >-HHAX 
LMAX-M-1 

2B  J-B 

DO  IB  I-l.LMAX 

IF  <H<  I*  1  )  CE  H(  I)>  GO  TO  IB 

J -J  ♦  I 

DUM-H( IM> 

H< 1*1  )-N<  I  ) 

H< I )-DUM 
IB  CONTINUE 

IF  (J  NED  GO  TO  2B 
I  *■ 

41  I-IM 

33  IF  (I  CT  LMAX )  GO  TO  3B 

IF  <H(  !♦!  )  HE  H<  I  )  )  GO  TO  4B 
DO  SB  J-I,LMAX 
SB  H<  J  >-H<  J*1  ) 

LMAX-LNAX-1 
GO  TO  39 
3B  IH-LMAXM 
RETURN 
END 

SUBROUTINE  SELECT  < X , Y, Z. KMAX , U , V , W , MAX, J , WW ) 

C 

C  THI 8  SUBROUTINE  SELECTS  DATAPOINTS  WITH  A  GIVEN  X  CO-ORDINATE 

C 

DIMENSION  X<  KMAX >, Y<  KMAX), ZC KNAX),V( MAX), W( MAX  ) 

I 

J-B 

EPS-ABSC  B  1-WW) 

IB  1-1*1 

IF  (I  CT  KMAX  )  GO  TO  2B 

IF  <ABS<  X<  I  >-U>  GT  EPS)  GO  TO  lfl 

J-J*l 

V<  J >-V<  I  ) 

W(  J  >-Z<  I  > 

CO  TO  IB 
2B  RETURN 
END 

SUBROUTINE  ORDER  (V,W,LMAX> 

C 

C  THIS  SUBROUTINE  0RDER8  THE  VALUES  OF  V  INTO  ASCENDING  ORDER 

C 

DIMENSION  V<  LMAX  ) , W( LMAX  ) 

IF  (LMAX  EM)  RETURN 
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LLMAX-LMAX-1 
28  J-B 

DO  IB  I-  I »  LLNAX 

IP  ( V( I* 1 ) . GE . V( I >  )  GO  TO  IB 

DUN*V( 1*1  ) 

V<I*l  )-V<I  ) 

V(I  )*OUH 
DUH -U< 1*1 > 

W<I*1  )-U< I  ) 

W< I )-DUH 
IB  CONTINUE 

IF  (J.NE.I)  GO  TO  2B 

RETURN 

END 

SUBROUTINE  CMECK2  < V » LMAK • HU  ) 

C 

C  THIS  SUBROUTINE  EHSURE8  THAT  IH  THE  SITUATION  WHERE  A  GRID  LIHE  JUST 
C  TOUCHES  A  CONTOUR  THE  TWO  DATA  POINTS  REPRE3EHTIHG  THE  POINTS  OF 

C  INTERSECTION  OF  THE  GRID  LINE  WITH  THE  CONTOUR  ARE  NOT  COINCIDENT 

C 

DIMENSION  Y<  LMAX  > 

IF  <LMAX.LT. 2)  RETURN 
EP8*B . B1*WW 
1-1 

IB  1-1*1 

IF  <!  CT  LMAX)  RETURN 
IF  <  Y< I ) . HE . Y< I - 1 >  >  CO  TO  IB 
IF  (I.GT.2)  YL-Y< I -2  ) 

IF  (I.LT.LMAX)  VH-Y< 1*1) 

IF  < I . EB . 2  )  YL-Y< 1-1  )-3  B*EPS 
IF  < I . EB .LMAX  )  YH-Y< I >*S . B*EPS 
88-AB8C YL-Y< 1-1  ))/3.B 
IF  < S8  . GT . EPS  )  SS-EPS 

Y<I -l  )-<  < YL-Y< 1-1 )  )/AB8<YL-Y<I-l))*SS>*Y< 1-1 > 

SS*ABB<YH-Y< I ))/9. I 
IF  <88.CT.EPS)  88«EP8 

Y< I  )-<  < YH-Y< I >>/ABS< YH-Y(  I )))*8S*Y< I ) 

CO  TO  IB 
END 

SUBROUTINE  AVTAN  < A. Y ,1 , N , YY « 22 , JOIN , H , I H , FRAC TH , I  ROUTE / 2E ) 

THIS  SUBROUTINE  CALCULATES  EQUISPACED  DATA  POINTS  ALONC  A  CRID  LINE 

DIMENSION  A<N).YY<  JOIN  ),ZZ<  JOIN).  Y<H),Z(N  >,H<  IH) 

El-Z< 1  ) 

E2-Z< N ) 

IF  < I ROUTE . NC . 1 )  GO  TO  3 

CI.W  ' 


3  IF  <H.EB  1)  60  "0  7B 

A<  1  >-<Z<2)-Z<  1)  V<  Y<2)-Y<  1  )) 
A<N)-<Z<N)-Z<N-i  >)/<  Y<N)-Y<M-I  )  ) 

IF  <N . CO .2 )  GO  TO  13 

INT-N-1 

DO  IB  1-2. INT 

A1-<Z< I >-Z< I-l>>/< Y< 1  )-Y< 1-1  )) 
A2-<  Z< I*l)-Z< !))/< Y< 1*1 )-Y< I >> 

IB  A<I  )-<Al*A2)/2  8 
13  J*B 


C 

C 

c 
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1-1 

20  J -J  ♦  1 

IF  <d  GT  JOIN)  RETURN 
IF  (  YY< J  )-Y< I  ))  25/30/30 
25  2 Z<  J  ) * E 1 
GO  TO  20 
35  J-J*l 

IT  < J  GT  JOIN  >  RETURN 
30  IF  <YY<J)-Y<I>)  20,40/45 
40  ZZ<  J  )-Z(  I  ) 

CO  TO  35 

45  IF  < YY< J >-Y< 1*1  )  )  55.60/65 
65  1-1*1 

IF  < 1*1  . GT . N)  GO  TO  50 
GO  TO  30 
60  ZZ<  J)-Z<  1*1  ) 

GO  TO  35 

55  CALL  AVTAN2  ( Y< I  ) , Y< I ♦ 1  ) , Z < I  > , Z< I ♦ 1  )  , A< I  )  . A( I ♦ 1  ) , Y Y< J  ) , Z Z< 0  ) 

• , H, IH/FRACTN) 

CO  TO  35 
70  J-l 

50  DO  75  L-J/JDIH 
75  ZZ<  L  )-E2 
RETURN 
END 

SUBROUTINE  AVTAN2  < Y 1 , Y 2, Z 1 , Z2 . A1 , A2 , Y Y , Z Z . H , I H , FR ACT N > 

THIS  SUBROUTINE  IS  THE  AVERAGE  TANGENTS  ALGORITHM  FOR  CONTOUR  CUT  DATA 
IT  ENSURES  THAT  INTERPOLATED  VALUES  REMAIN  UITHIN  THE  CONTOUR  INTERVAL 
UNDER  CONSIDERATION  OR  FOR  TURNING  POINTS  DOES  NOT  DEVIATE  FROM  THE 
INITIAL  CONTOUR  DATA  BY  MORE  THAN  ONE  CONTOUR  INTERVAL 

DIMENSION  H< I H ) 

DUN  1* A 1 
DUM2-A2 
DUM3-Z 1 
DUM4-Z2 
DUM5-YY 

IF  (Z2.CE.Z1)  GO  TO  5 
DUM-Z1 
Z1-Z2 
Z2-DUM 
DUM-A 1 
A1--A2 
A2--DUM 
Y Y-Y2+Y1 -YY 
5  DY-Y2-Y1 
DZ-Z2-Z1 

A*< Al*0Y*A2*DY-2  .  0*OZ  )/<  DY**3> 

B»<  3. 0*DZ-2  0*A1*DY-A2*OY  V<  DY**2  > 

C  -A  1 
Y-YY-Y1 

ZZ-A*Y**3*0*Y**2*C*Y*Z1 
IF  <Z2  HE  21 >  CO  TO  IB 
IF  ( A2  .  EQ  A  1 )  CO  TO  70 
10  I  TEST -0 

TEST-B**2-3  0*A*C 

IF  <  A8S<  A*DY* 150  0  >  GT  A0S<B>>  GO  TO  15 
YS1— C/(  2  0*0  ) 

IF  ( YS1  LT  0  0  OR  YS1  GT  DY)  GO  TO  70 
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ITEST-1 
GO  TO  35 

15  IF  (TEST  )  7  0/ 28/ 25 

2B  YSl m-B/l 3 . B*A  ) 

IF  (YSl  LT.B.B.OR  YSl  GT.DY)  GO  TO  7 B 

ITEST-1 

GO  TO  35 

25  YSl-< -8-TEST**0 . 5  i/<  3  .  P*A  ) 

IF  < YSl  . LT  B  B  OR. YSl  .GT. DY>  GO  TO  3B 
ITEST-1 

3B  YS2»< -B*TEST**0  5>/(3.0*A> 

IF  (YS2.LT.8.B.0R. YS2.GT.0Y)  GO  TO  35 

IF  (ITEST  EC  B)  YS1-YS2 

ITEST-ITEST+1 

IF  < ITEST  HE  2)  CO  TO  35 

IF  ( YSl  LT  YS2)  CO  TO  36 

DtlM-YS  1 

YS1-YS2 

VS2-0UM 

36  IF  ( A 1  HE  B.R)  GO  TO  3? 

YS1-YS2 

ITEST-1 

IF  ( A2  HE  8 .0  )  GO  TO  35 

ITEST-B 

GO  TO  35 

3?  IF  vA2  HE  1.8)  GO  TO  35 
I  TEST-1 

35  IF  ( ITEST  EC.  0)  GO  TO  70 
IF  < ITEST  EC. 2)  GO  TO  55 
1  SI GH- ♦ 1 

IF  (A1-A2  LT.B.B)  ISIGH--1 
DO  4B  1-1.  IH 

4 B  IF  <Z1  EO  H< I  ))  GO  TO  45 

45  IF  ( ISICN  Ell)  I-I+l 
COHT-HC I  )-H( I  - 1  ) 
ZT-A*YS1**3*8*YS1**2*C*YS1 
IF  (ABS(ZT).LT.CONT)  GO  TO  ?B 
DZ-FLOAT( IS1GH)*FRACTN*C0NT 
IF  ( Y  CT  YSl  )  GO  TO  5fl 
A -( Al* YSl -2 .0*02  >/ YS 1 •• 3 
B»(  3  0*02-2  0*A1*YS1  >/YSl**2 
22-A*Y**3*B*Y**2*C*Y+Zl 
IF  ( ABS( ZZ-Z1 >  GT  ABS(OZ)  )  ZZ-Z1*DZ 
GO  TO  7B 

5B  DY-DY-YS1 
Y-Y-YS1 

A-( A2  *  0 Y  *2 . 8*02  )/0Y**3 
B -< -3  B*D2-A2*0Y  VDY**2 
ZZ-A*Y**3*B*Y**2*Z1*0Z 
IF  (ABS<ZZ-Zi  )  GT.ABS(OZ>)  22*21*02 
GO  TO  78 

55  ZT1-A*YS1**3*8*YS1**2*C*YS1*Z1 
ZT2»A*YS2**3*B*YS2**2*C*YS2*ZI 
IF  (2T1.LT  22. AMD.ZT2.GT  Z 1 >  CO  TO  7B 
IF  (Y  CT  YSl  )  CO  TO  6B 
IF  (2TI.LT  22)  GO  TO  70 
DZ-FR ACTH*( Z2-21  ) 

A»<  Al *YSl-2 . B*0 2  )/YSl**3 
B-( 3  8*02-2  0*A1*YS1 >/YSl**2 
ZZ*A*Y**3*8*Y**2*A1*Y*Z1 


noon 
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IF  (22  GT  Zl  +  DZ  )  ZZ-Z1+DZ 
GO  TO  ?fl 

6 0  IF  (Y  CT  YS2)  CO  TO  65 

IF  CZT1  GE  Z2  )  ZT1-Z1 +FRACTN*( Z2-Z1  > 

IF  (  Z  T  2  LE  ZI  )  ZT2»Z2-FRACTN*< Z2-Z1 > 

DZ-2T2-ZT1 

DY-YS2-YS1 

Y*Y-YS1 

A*-2  B*DZ/D Y**3 
8=3  B*DZ/DY**2 
ZZ-A*Y**3+B*Y**2+ZT1 
IF  (ZZ  CT  ZT1 >  ZZ=2T1 
IF  (ZZ  LT  ZT2  )  Z2-2T2 
CO  TO  7B 

65  IF  <2T2. CT .21)  CO  TO  7fl 
0Z“FR ACT  N*< 22-21  ) 

DY-DY-YS2 

Y-Y-YS2 

A-( A2  *  DY -2 B*DZ  >/DY**3 
B»( 3  0*DZ-A2*OY  >/0Y**2 
2Z*A+Y**3+0*Y**2+Z2-DZ 
IF  (ZZ  LT  Z2-OZ)  22*Z2-DZ 
7B  A1-OUH1 
A2-DUH2 
21-DUM3 
Z2-DUM4 
YY*DUM5 
RETURN 
END 

SUBROUTINE  «iICHT  ( Y , Y Y , U , N , JD I M , Z . A > 

THIS  SUBROUTINE  GENERATES  WEIGHTING  VALUES  ASSOCIATED  WITH  THE 
INTERPOLATED  DATA  POINTS  ALONG  A  GRID  LINE 

DIMENSION  YY(JDIN),W(JDIfl),Y(N>,Z(N>,A(N) 

J  =  1 

IF  (N  EQ  1  )  CO  TO  78 
J  =B 
1*1 

IB  J=J*1 

IF  ( J  GT  JDIM  )  RETURN 
IF  (YY(J  )- Y ( I  ))  2B , 38 , 3B 
28  W( J  >  =  B  B 
CO  TO  IB 
4 B  J*J*1 

IF  ( J  CT  JOIN )  RETURN 
3B  IF  ( YY< J >-Y( I  +  l  ) )  38,50,68 
6B  1-1*1 

IF  ( I+l  CT  N)  GO  TO  7B 
GO  TO  3B 

5B  W<  J  )»  ABS< 1  B/( Y(  1*1  )-Y(  I  )  )  ) 

IF  ( Z(  I )  EQ  Z( I *1 )  AND  A(  I  )  EQ  B  0  AND  A(  I ♦ 1 >  EQ  8  8 )  W(J>=-1  8 
GO  TO  4B 

7B  DO  88  L-J.JDIN 

88  W(L  )-B  0 
RETURN 
END 
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SUBROUTINE  COORDS  C F , G, H , X , Y - Z ,  I H AX , KH AX , WF , X ,  I E X I T > 

THIS  SUBROUTINE  EVALUATES  THE  CO-ORDINATES  AT  WHICH  THE  GRID  LINES 
INTERSECT  A  CONTOUR.  LINEAR  INTERPOLATION  IS  EHPLOYED 

D 1HENS ION  F( I  MAX  >, C< I  HAH), X(KHAX), V( XRAX  >,Z(XHAX  > 

EPS»A8S<  B . 0 1*MF ) 

1-3 

I  FI  -NF  IX<  F<  I-D/UF) 

IF  (  AB$(  FLO AT < IF  1 >*MF-F( 1-1 ) >. CT . EPS  )  GO  TO  5 

X-X  ♦  1 

X(  X  )-F(  I  ) 

Y  (  X  >«G(  I  ) 

Z(X  >-H 
5  I-I-l 
88  I » I ♦ 1 

IF  (  I  .  CT  IHAX-1  )  RETURN 
IF1-HFIX(F(  I-D/UF) 

I F2-NF I X(  F (  D/UF  > 

FF-(F(  1  >-F<  I)  )■*<  F(  I  )-F(  1-1  )> 

NF-IF2-IF1 
IF  <  NF  )  4B, 58,68 
68  DO  28  N- 1 , NF 
X  -X  ♦  l 

X<  X  >-FLOAT  <  IFl)*WF+'FLOAT<N)*UF 
IF  (ABS( X( X  )-F( I  )>  CT  EPS )  CO  TO  18 
Y(X  )-G(  I  ) 

Z(X  )-H 

IF  (FF)  118,88,88 

18  Y(K >-( (X<X  >-F< 1-1 >>•( G< I )-C< 1-1 >>A FC I >-F< 1-1 >  ))*C<  1-1  ) 

28  Z(  X  )-H 
GO  TO  88 
48  NF--NF 

IF  (  ABS(  FLOAT( IF2)*8F-F( I ) >. LT . EPS)  NF-NF+1 
DO  38  N-l.HF 

XX-FLOATC IF1 )*BF-FLOAT( N-l  )*UF 

IF  < ABS( XX-F( 1-1  )>  LT  EPS  )  GO  TO  3B 

X-X*l 

X(X  )«XX 

IF  CABS< X< X  )-F< I  ))  CT  EPS )  GO  TO  78 
Y<  X  >-G( I ) 

Z(X  )-H 

IF  (FF )  118,88,88 

78  Y(X  )«(  <X(X )-F( 1-1 ) >*<  G< I )-G<  1-1  >>/( F<  I  )-F( 1-1  )  )  >+G(  1-1  ) 

3B  Z(X)«H 
GO  TO  88 

58  XX-FLOAT( IF1 >*MF 

IF  <ABS(XX-F( I))  GT.EPS)  GO  TO  88 
X-X  +  l 

IF  (X  CT . XH AX  )  GO  TO  128 
X(X  )-F( I ) 

Y(X  >»G( I ) 

Z<X  )«H 

IF  (FF )  118,88, 88 
118  K-K+l 

X(X  )«X(X-1  ) 

Y(X  )-Y(X-l ) 

Z(X  >«Z(X-i  ) 

GO  TO  8B 
ENO 
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El 

C 

C  •••••  PROGRAM  DRIVER1  ••••• 

C 

C  THIS  PROGRAM  PRODUCES  CONTOURS  OF  «  2-0  CAUS8 ! AM  >  CENTRE  <■.■> 

C 

DIMENSION  N<13),Y<  13) 

I  TYPE  -  8 
INC-12 

Pl-3  1 41 392494 
H  MAX- 2  40  0 
HPBU-18  8 
CONT-23  8 

CONST  -AL0C<  2 . 8  )/'MPBM--2 
DO  18  J-1,8 
F-CONT-FCOAT<  IB- J  ) 

R-(  <  ALOGC  MHAXXF  )  l/COHST  )Ml.l 
DTHETA-2  8-PI /FLOAT! INC  > 

DO  28  I-l. INC 
THETA-DTHETA-FLOAT( t-1  ) 

N (  I  )-R-C08(  THETA  ) 

28  Y < I  )-R-S INC  THETA  ) 

INUM-INC-l 
X< INUH  >•*<  1  > 

Y< INUH  )- Y<  1  ) 

WRITE  <8.81  )  F. INUN. I  TYPE 

81  FORMAT  <  F8  1. IS. 12  > 

WRITE  (8.82)  (N<  I  ),Y<  1  >.  I-l.  INUH) 

82  FORMAT  <2F8  3) 

18  CONTINUE 

STOP 

END 


E2 


C  •••••  PROCRAM  DRIVER2  ••• •• 

C 

C  THIS  PROGRAM  EVALUATES  A  2-0  GAUSSIAN  OVER  A  38-38  CRID 

C 

DIMENSION  2(38) 


PI-3  141382834 
HMAX-248  8 
HPBW-18  ■ 

C0NST-AL0G<2  8  >/HPBW--2 

IDIM-38 

JDIM-38 

WRITE  <8.88)  IDIM.  JOIN 
88  FORMAT  < 13, IN. 13  ) 

00  18  J-l.JDIM 
Y  — 28  0*FLOAT<  J-l  ) 

DO  28  I-l, IDIM 
N  —  28  B-FIOATCI-I) 
RR-K--2-Y--2 

2<I)-24B  B/EXP<  C0N8T-RR  ) 

28  IF  (2(1)  IT  23  8  )  2<  I  )-23  I 
18  WRITE  (8.81)  <2(  1), I-l.  38) 
81  FORMAT  <  1 2 F 8  1  ) 

STOP 

END 
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Appendix  F 
PROGRAM  EXPAND 


>1 


F . 1  Function 

To  expand  (or  contract)  the  number  of  points  in  a  grid  in  either  or  both  the  x 
and  y  directions. 

F. 2  Dimensions  of  arrays 

number  of  points  in  a  row  or  column,  whichever  is  greatest  (ID1M) 

number  of  points  in  a  column  after  expansion  (JJDIM) 

TDIM,  JJDIM. 

F. 3  Program  units 

Subroutine  AVTAN  :  Calculates  equispaced  data  points  along  a  grid  line.  (This 
subroutine  is  very  similar  to  that  described  in  Appendix  B.  The  main  difference  is  that 
interpolated  values  outside  the  dataset  are  only  set  to  the  value  of  the  nearest  point  in 
the  dataset,  no  other  option  being  available.) 


X 
A 

Z 

XX  ) 
ZZ  f 

AA 


Subroutine  AVTAN 2  :  Average  tangents  algorithm  for  uniformly  spaced  data. 

In  addition  the  following  FORTRAN  functions  are  called: 


IF1X  :  integer  value  nearer  to  zero 
FLOAT  :  converts  integer  to  floating  point. 


F.4  Input/output 

Input  via  FORTRAN  channel  I 
Input  via  FORTRAN  channel  5 
Output  via  FORTRAN  channel  6 

F. 5  Listing  of  program  EXPAND 


expansion  factors  for  x  and  y  directions 
the  input  grid  (formatted  as  output  by  MAPGRID) 
expanded  grid  (format  similar  to  that  of  input 
grid)  . 


1  10 
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C  PACE  Mil 

C 

C  PROGRAM  E  X  PA  HO  ••••• 

C 

c  this  procrah  expands  or  contracts  the  number  of  points  in  a  grid 


c 

COMMON  X(237).XX(  H3S  >.  A(  237  >,  2(  237  )  .  ZZ(  1 193 >  .  RA(  297 .  1033  > 

COMMON  AT<  237.297) 

REAO  (1.11)  KM.  YH 
1  1  FORMAT  <  2FB  3  ) 

REAO  (3.31)  IDIM.JDIH 

31  FORMAT  ( 13. IX. 13) 

DO  10  J-l.JDIM 

10  READ  (3.32)  (  AT  (  I ,  J  ) ,  !•  1 .  1 01 M  > 

32  FORMAT  (  I2FE  1  ) 

IPTS-IFIXC  XN-FLOAT( IOlM-1 >  >♦! 

JPTS-1FIX<YN-FL0AT(J0IN-1>>»1 
WRITE  ( 6 . E 1  )  IPT8.JPT8 
El  FORMAT  ( 14. IX. 14  > 

SCALE-FLOAT! J DIN-1  V FLOAT! JPTS- I  ) 

DO  20  J-l, JPTS 
20  XX(  J  )-FLOAT(  J-l  )*8CALE 
DO  30  J-l. JOIN 
30  X( J  )-FLOAT( J-l  ) 

DO  40  I-l.  I  D  I  M 
DO  30  J-l. JOIN 
30  2(  J  )•  A  T(  I.  J  > 

CALL  AVTAH  (A, X.  2,  JOIN. XX. 2Z, JPTS) 

DO  40  J-l, JPTS 
40  AA(  I,  J  >-2Z(  J  ) 

BCALE -FLOAT ( IDIH-1  VFLOAT! IPT8-I  ) 

DO  E8  I-l. IPTS 
E0  XX< I )-FLOAT( I-l  )*BCALE 
DO  78  I-l. ID1H 
70  XI  I  J-FLOAT (  I-l  ) 

DO  80  J-l, JPTS 
DO  90  I-l. IDIH 
90  Z(I  )-  A  A<  1,  J  > 

CALL  AVTAN  ( A , X , Z. ID  I M, XX, 22 . IPTS > 

80  WRITE  (  E  .  32  >  (ZZ(1  ).  I-l. IPTS) 

STOP 

END 

SUBROUTINE  AVTAH  ( A, V . Z . H . TV . Z Z , JO IIU 
C 

C  THIS  SUBROUTINE  CALCULATES  E0U I  SPACED  ORTA  POINTS  ALONG  A  CRIB  LINE 

C 

DIMENSION  A(N).YY(JOIM).ZZ(JOIN).T<N).Z<N) 

£1-2(1) 

E2*Z<  N  ) 

IF  (M  EO  1  )  CO  TO  70 

A<  I  >-(Z(2)-Z<  1  >)✓(  Y<  2  )-Y(  1  >> 

A(N  >•(  ZC  N)-Z(  N-l  >>/!  V(N  >-Y(N-l  )  ) 

IF  (H  LT  3)  GO  TO  15 

INT-N-1 

00  10  1-2. INT 

AI-(Z(  I  )-Z(  I-l  >  )/(  Y<  I  »-VC  I-I  >> 

A2-(Z(  1*1  >-Z(  I  )  >7!  Y(  I-l  )-Y(I  >) 

10  A<  I  )•(  A1  «A2  >72  0 
13  J-B 
1-1 

20  J - J ♦  I 
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C  PACE  BBB2 


IF  (J  CT  JOIN)  BE  TURN 
IF  <YY(J)-Y<I>)  29,31,31 
29  Z2(  J  )-EI 
CO  TO  21 
39  J-J*l 

IF  (J  CT  JOIN)  RETURN 
3>  IF  <YY(J)-Y<I>)  29,41.49 
41  ZZ<  J>-Z(  I) 

CO  TO  39 

49  IF  < YYC J >-YC 1*1  )  )  99.<l.tS 
49  1-1*1 

IF  <1*1  CT  N>  CO  TO  SI 
CO  TO  31 
Cl  2Z<  J)-ZC  1*1  ) 

CO  TO  33 

99  CALL  AVTAN2  <  AM  *1  ).  A<  I  ),  YY<  J  ) .  22<  J  ) .  Y<  I  *1  >.Y  <  I  )  .  2<  I*  1  > .  Z<  I  >  > 

CO  TO  33 
71  J-l 

31  DO  79  L-J. JOIN 
79  22<  L  >-E2 
RETURN 
END 

SUBROUTINE  AVTAN2  < A2 . A 1 , Y Y . ZZ . Y2 . Y I . 22. 2 1 > 

C 

C  THIS  SUBROUTINE  II  THE  AVERAGE  TANGENTS  ALGORITHN  FOR  EVENLY  SPACED  DATA. 

C 

OY-Y2-YI 

02-22-21 

A-< Al*0V*A2«0Y-2  1*02  )7<0Y««J> 

B-< 3 . 1*02-2  B-Al »0Y-A2*DY )/( 0Y»«2  ) 

C-Al 

Y-YY-YI 

ZZ-A-Y  — 3*S-Y  — 2*C-Y*2l 

RETURN 

ENO 
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Appendix  C 
PROGRAM  SECTION 

G. I  Function 

To  produce  a  cross-section  across  a  grid  of  points. 

G. 2  Dimensions  of  arrays 

X  ) 

^  I  maximum  number  of  points  required  in  a  cross-section 
IW  ] 

ZT  4 

ZZ  number  of  points  in  a  row  in  the  grid,  4. 

G . 3  Program  units 

Subroutine  AVTAN  :  average  tangents  algorithm  for  uniformly  spaced  data. 

In  addition  the  following  standard  FORTRAN  functions  are  called: 

FLOAT  :  converts  integer  to  floating  point 
IFIX  :  integer  value  nearer  to  zero. 

G.4  Input/output 

Input  via  FORTRAN  channel  1  :  Co-ordinates  of  top  left  corner  of  grid  and  grid 

spacing. 

Co-ordinates  of  the  two  end  points  of  the  required 
cross-section. 

Number  of  points  required  in  the  cross-section. 
Input  via  FORTRAN  channel  5  :  The  input  grid  (formatted  as  output  by  MAPGRID) . 

Output  via  FORTRAN  channel  6  :  Number  of  points  in  the  cross-section  and  the 

x  and  y  co-ordinates  and  the  height  of  the 
points  in  the  cross-section. 

G. 5  Listing  of  program  SECTION 


o  r>  n  o  o  o  r>  o  o  ooor>  non  oor>p»  o  o  o 


Appendix  G 


c 


PAGE  0001 


•«*•»  PROGRAM  SECTION 

01  HENS  ION  ZZ(  1833/ 4 ) »  X( 2B8B  ) . Y( 2IBB).Z(2BBB  >. ZT( 4  ) » IUC2BBB) 

INPUT  CO-ORDINATES  OF  TNE  TOP  LEFT  NANO  CORNER  OF  GRID  AND  THE  GRID 
SPACING 

READ  <1.11)  XI. YB.W 

INPUT  THE  CO-ORDINATES  OF  THE  TWO  ENDS  OF  THE  REQUIRED  CROSS-SECTION 
READ  < 1. 12)  X1.V1.X2.Y2 


INPUT  NUMBER  OF  POINTS  REQUIRED  IN  THE  CROSS-SECTION.  IF  THE  NUMBER  IS 
LES8  THAN  2  THE  PROCRAH  CHOOSES  A  NUMBER  BASED  ON  THE  GRID  SPACING 

READ  ( 1. 13)  IPTS 

CHANCE  DIRECTION  OF  CUT  IF  Yl.LT.Y2 
IDIRN-1 

IF  < Y1  .GE.  Y2)  GO  TO  3 
DUN-X1 
X1-X2 
X2-DUM 
DUM-Y1 
Y1-Y2 
Y2-0UM 
IDIRH--1 
3  CONTINUE 

CALCULATE  P0INT8  AT  WHICH  INTERPOLATED  VALUES  ARE  REQUIRED 

IF  < IPTS  GE .2)  GO  TO  18 
D-<<X2-X1)**2+<Y2-Y1  )**2)**B  3 
IPTS-IFIX<D/W)«-1 
11  3X«<X2-X1 )/FLOAT< IPTS-1 > 

SY«<Y2-Y1VFL0AT<IPTS-1  ) 

DO  13  1-1. IPTS 
X< I  )»X1*FL0AT< 1-1  )*SX 
13  Y< I )«Y1*FL0AT< 1-1  )*8Y 

PERFORM  INTERPOLATION 

READ  (3.31)  IDIH.JDIH 
DO  28  K-1.4 

28  READ  (3.32)  < ZZ< I. K ). I- 1 . IDIH ) 

J«2 

Y8“YB*FL0AT< I  DIM-1  )*U 
DO  23  L* 1 . IPTS 
IX-IFIX<<X(L)-XB)/W)*1 
I Y“ IF  I X<  <  Y8-Y<  L  )  )/W)M 

XREL-(<X(L  )-XI)/W)-FLOAT<  I  FI  X<  <  X<  L  )-XB  VW  >  ) 

YREL«<  <Y8-Y<L  ))/W> -FLOAT!  IFIX<  <  YB-Y<  L)  VW  )> 

IF  < IX  LT  2  )  CO  TO  33 
IF  (IX  CT. IDIM-2 )  GO  TO  33 
IF  (IY.LT  2)  CO  TO  33 
IF  <!Y  CT  JDIM-2)  CO  TO  33 
38  IF  < I Y  EQ  J  )  GO  TO  48 


Appendix  C 


-•  i 


MCI  hi: 


c 

c 

c 


c 

c 

c 


00  31  K*l,  3 
00  33  1*1.  tom 
33  2  2<  I .  K  )- 22  <  I.K*I  > 

J.J»| 

READ  (3.32)  <  Z2<  1.4),  1*1,  IOIN) 

CO  TO  3fl 
4  B  00  43  X*  1 . 4 

21*22<  I X -2  *X , 1  ) 

Z2*Z2< IX-2*X, 2) 

23  *22<  I X  -  2  ♦  X  ,  3) 

Z4*22<  1X-2*X, 4) 

43  CALL  AVTAN  <  2  1 . 22, 23, 24  , YREl ,ZT<* >> 

CALL  AVTAN  (2T(  I  >,  2T<  2),ZT<  3  ),  2T<  4). KIEL.  ZT  I  > 

00  3i  K*1 , 4 
2  1  *  22<  IX-I ,K> 

22*221  IX, K  > 

2  3*  22(  I  X ♦ 1  ,  X  ) 

Z4*Z2<  I  X *2 , X  ) 

SB  CALL  AVTAN  (21 . 22. 23, 24 ,XBEL ,2T(X )) 

Call  AVTAN  <  ZT<  1  >,  ZT<  2 ) ,  ZT<  3  >.  2T<  4  > ,  YREL ,  ZT2  > 

2 (  L  >■(  ZT l»2T2  >72  ■ 
t  W(  L  )*1 
CO  TO  23 
33  2(L  >*B  B 
t  W(  L  >*B 
23  CONTINUE 

OUTPUT  DATA  (DEPENDS  UPON  MHETHER  THE  DIRECTION  OF  TNE  CUT  MAS  REVERSED) 
L*B 

AB  L  *L  ♦  1 

IF  ( IM(L  )  EO  fl>  CO  TO  4B 
CPAD* Z( L ♦ 1  )-2(L  ) 

LNAX-L-I 

DO  AS  LL-1 . LNAX 

A3  Z(LL)*2(L)-FL0AT(L-LL)*CRAD 
7B  L  *L ♦ 1 

IF  ( IM(L  >  EO  I )  GO  TO  7B 
GRAD*  Z( L - l  )-2(L-2) 

DO  73  LL*L ,  IPT8 

73  2(  LL  >*Z(  L-  1  )-CRAD*FLORT(Ll-L*l  ) 

WRITE  (A. 41)  IPTS 
IF  < IDIRN  EO  I)  CO  TO  S3 
DO  Bfl  1-1. IPTS 
I I*IPTS*1- I 

SB  WRITE  (  A  ,  A  2  )  X(  I  I  ) ,  V  (  1 1  >,  2  (  I  I  > 

STOP 

83  00  SB  1*1 . IPTS 

SB  WRITE  (A. 42)  X(  I  >,  V(  I  ),  Z(  I  ) 

STOP 

INPUT /OUTPUT  FORNATS 


1 1 

FORMAT 

( 3F6  2 ) 

12 

FORMAT 

(  4F6  2  > 

1  3 

FORMAT 

<  14  ) 

51 

FORMAT 

(  I  3  r  1 X  #  I  3  > 

52 

FORMAT 

(  1 2  F  6  1  > 

61 

FORMAT 

<  14  ) 

62 

FORMAT 

(  3F6  1  > 

END 

SUBROUTINE  AVTAN  ( 2 1 , 22 , 2 3 . 2 4 , V Y , 22 ) 

C 

C  THIS  S  IP  OUTINE  IS  THE  AVERACE  TANCENTS  AL C OR  I THN  FOR  UNITY 

c  spaceo  ; a r a  points 
c 

Al* <23-21 >72  B 
A2*<  24 -22  >72  B 
02*23-22 

A*(  At  *A2-2  B*D2  ) 

B-<  3  B*02-A2-2  8*A1  ) 

C*A  1 

22*A*YY**3«8*YV**2«C*YY*22 

RETURN 

END 
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Numbers  shown  represent  (grid  value  ;  25)  rounded  to  the  lower  integer 


Fig  4  Grid  produced  by  taking  cuts  in  y  direction 
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Fig  5 
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Numbeis  shown  represent  (quit  value  751  rounded  to  the  lower  mteqei 


Fig  5  Grid  produced  by  taking  cuts  in  x  direction 
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Fig  6  Magnitude  of  the  differences  between  the  program  and  the  analytic 
grids  using  method  1  section  4.4 
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Fig  7  Magnitude  of  the  differences  between  the  program  and  the  analytic 
grids  using  method  2  section  4.4 
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Fig  8  Magnitude  of  the  differences  between  the  program  and  the  analytic 
grids  using  method  3  section  4.4 
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Numlrers  shown  represent  (qriri  value  ?!>)  rounder!  to  rhe  lower  integer 


Fig  10  Final  grid  produced  by  combining  those  shown  in  Figs  4  end  5 
using  method  3  section  4.4 
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Numbers  shown  represent  (weighting  value  >20)  rounded  to  the  lower  integer 


Fig  13  Weighting  values  associated  with  each  point  of  the  final  grid 


Fig  14a-d  Operation  of  program  MAPGRID  on  a  contour  map  consisting  of  two  sets  of  contours 


a  Original  contour  map 


—  False  contour  0<height<1  - - 

b  Contour  map  with  false  contours  added 


Fig  15a&b  Introduction  of  false  contours 
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Fig  16  Geometry  of  the  interpolation  performed  by  program  SECTION 
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